John von Neumann Institute for Computing 



Ab Initio Molecular Dynamics: Theory and 
Implementation 

Dominik Marx and Jurg Hutter 



published in 

Modern Methods and Algorithms of Quantum Chemistry, 

Proceedings, Second Edition, J. Grotendorst (Ed.), 

John von Neumann Institute for Computing, Julich, 

NIC Series, Vol. 3, ISBN 3-00-005834-6, pp. 329-477, 2000. 



© 2000 by John von Neumann Institute for Computing 
Permission to make digital or hard copies of portions of this work for 
personal or classroom use is granted provided that the copies are not 
made or distributed for profit or commercial advantage and that copies 
bear this notice and the full citation on the first page. To copy otherwise 
requires prior specific permission by the publisher mentioned above. 



http://www.fz-juelich.de/nic-series/ 



AB INITIO MOLECULAR DYNAMICS: 
THEORY AND IMPLEMENTATION 



DOMINIK MARX 

Lehrstuhl fur Theoretische Chemie, Ruhr-Universitat Bochum 
Universitatsstrasse 150, 44780 Bochum 
Germany 

E-mail: dominik. marxQtheochem. ruhr-uni-bochum. de 

JURG HUTTER 

Organisch-chemisches Institut, Universitdt Zurich 
Winterthurerstrasse 190, 8057 Zurich 
Switzerland 
E-mail: hutter@oci.unizh.ch 

The rapidly growing field of ab initio molecular dynamics is reviewed in the spirit 
of a series of lectures given at the Wintcrschool 2000 at the John von Neumann 
Institute for Computing, Jiilich. Several such molecular dynamics schemes are 
compared which arise from following various approximations to the fully coupled 
Schrodingcr equation for electrons and nuclei. Special focus is given to the Car— 
Parrinello method with discussion of both strengths and weaknesses in addition 
to its range of applicability. To shed light upon why the Car-Parrinello approach 
works several alternate perspectives of the underlying ideas are presented. The 
implementation of ab initio molecular dynamics within the framework of plane 
wave-pscudopotcntial density functional theory is given in detail, including diag- 
onalization and minimization techniques as required for the Born-Oppcnhcimcr 
variant. Efficient algorithms for the most important computational kernel routines 
are presented. The adaptation of these routines to distributed memory parallel 
computers is discussed using the implementation within the computer code CPMD 
as an example. Several advanced techniques from the field of molecular dynam- 
ics, (constant temperature dynamics, constant pressure dynamics) and electronic 
structure theory (free energy functional, excited states) are introduced. The com- 
bination of the path integral method with ab initio molecular dynamics is presented 
in detail, showing its limitations and possible extensions. Finally, a wide range of 
applications from materials science to biochemistry is listed, which shows the enor- 
mous potential of ab initio molecular dynamics for both explaining and predicting 
properties of molecules and materials on an atomic scale. 



1 Setting the Stage: Why Ab Initio Molecular Dynamics ? 

Classical molecular dynamics using "predefined potentials", either based on em- 
pirical data or on independent electronic structure calculations, is well estab- 
lished as a powerful tool to investigate many-body condensed matter systems. 
The broadness, diversity, and level of sophistication of this technique is docu- 
mented in several monographs as well as proceedings of conferences and scientific 
schools 12 > 135 > 270 > 217 > 69 > 59 > 177 . At the very heart of any molecular dynamics scheme 
is the question of how to describe - that is in practice how to approximate - the 
interatomic interactions. The traditional route followed in molecular dynamics is to 
determine these potentials in advance. Typically, the full interaction is broken up 
into two-body, three-body and many-body contributions, long-range and short- 
range terms etc., which have to be represented by suitable functional forms, see 
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Sect. 2 of Ref. 253 for a detailed account. After decades of intense research, very 
elaborate interaction models including the non-trivial aspect to represent them 
analytically were devised 253 . 539 < 584 . 

Despite overwhelming success - which will however not be praised in this re- 
view - the need to devise a "fixed model potential" implies serious drawbacks, see 
the introduction sections of several earlier reviews 513 > 472 for a more complete di- 
gression on these aspects. Among the most delicate ones are systems where (i) 
many different atom or molecule types give rise to a myriad of different interatomic 
interactions that have to be parameterized and / or (ii) the electronic structure 
and thus the bonding pattern changes qualitatively in the course of the simulation. 
These systems can be called "chemically complex" . 

The reign of traditional molecular dynamics and electronic structure methods 
was greatly extended by the family of techniques that is called here " ab initio 
molecular dynamics" . Other names that are currently in use are for instance Car- 
Parrincllo, Hcllmann-Feynman, first principles, quantum chemical, on-the-fly, di- 
rect, potential-free, quantum, etc. molecular dynamics. The basic idea underlying 
every ab initio molecular dynamics method is to compute the forces acting on the 
nuclei from electronic structure calculations that are performed "on-the-fly" as the 
molecular dynamics trajectory is generated. In this way, the electronic variables are 
not integrated out beforehand, but are considered as active degrees of freedom. This 
implies that, given a suitable approximate solution of the many electron problem, 
also "chemically complex" systems can be handled by molecular dynamics. But 
this also implies that the approximation is shifted from the level of selecting the 
model potential to the level of selecting a particular approximation for solving the 
Schrodinger equation. 

Applications of ab initio molecular dynamics are particularly widespread in ma- 
terials science and chemistry, where the aforementioned difficulties (i) and (ii) are 
particularly severe. A collection of problems that were already tackled by ab initio 
molecular dynamics including the pertinent references can be found in Sect. 5. The 
power of this novel technique lead to an explosion of the activity in this field in terms 
of the number of published papers. The locus can be located in the late-eighties, 
see the squares in Fig. 1 that can be interpreted as a measure of the activity in 
the area of ab initio molecular dynamics. As a matter of fact the time evolution of 
the number of citations of a particular paper, the one by Car and Parrinello from 
1985 entitled "Unified Approach for Molecular Dynamics and Density-Functional 
Theory" 108 , parallels the trend in the entire field, see the circles in Fig. 1. Thus, 
the resonance that the Car and Parrinello paper evoked and the popularity of the 
entire field go hand in hand in the last decade. Incidentally, the 1985 paper by Car 
and Parrinello is the last one included in the section "Trends and Prospects" in 
the reprint collection of "key papers" from the field of atomistic computer simula- 
tions 135 . That the entire field of ab initio molecular dynamics has grown mature 
is also evidenced by a separate PACS classification number (71.15.Pd "Electronic 
Structure: Molecular dynamics calculations (Car-Parrinello) and other numerical 
simulations") that was introduced in 1996 into the Physics and Astronomy Classi- 
fication Scheme 486 . 

Despite its obvious advantages, it is evident that a price has to be payed for 
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Figure 1. Publication and citation analysis. Squares: number of publications which appeared 
up to the year n that contain the keyword u ab initio molecular dynamics" (or synonyma such 
as "first principles MD", "Car-Parrincllo simulations" etc.) in title, abstract or keyword list. 
Circles: number of publications which appeared up to the year n that cite the 1985 paper by 
Car and Parrincllo 08 (including misspellings of the bibliographic reference). Self-citations and 
self— papers are excluded, i.e. citations of Ref. 108 in their own papers and papers coauthored by 
R. Car and / or M. Parrinello arc not considered in the respective statistics. The analysis is based 
on the CAPLUS ("Chemical Abstracts Plus"), INSPEC ("Physics Abstracts"), and SCI ("Science 
Citation Index" ) data bases at STN International. Updated statistics from Ref. 405 . 



putting molecular dynamics on ab initio grounds: the correlation lengths and re- 
laxation times that are accessible are much smaller than what is affordable via 
standard molecular dynamics. Another appealing feature of standard molecular 
dynamics is less evident, namely the "experimental aspect of playing with the po- 
tential" . Thus, tracing back the properties of a given system to a simple physical 
picture or mechanism is much harder in ab initio molecular dynamics. The bright 
side is that new phenomena, which were not forseen before starting the simulation, 
can simply happen if necessary. This gives ab initio molecular dynamics a truly 
predictive power. 

Ab initio molecular dynamics can also be viewed from another corner, namely 
from the field of classical trajectory calculations 649 > 541 . In this approach, which 
has its origin in gas phase molecular dynamics, a global potential energy surface 
is constructed in a first step either empirically or based on electronic structure 
calculations. In a second step, the dynamical evolution of the nuclei is generated 
by using classical mechanics, quantum mechanics or semi / quasiclassical approx- 
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imations of various sorts. In the case of using classical mechanics to describe the 
dynamics - the focus of the present overview - the limiting step for large systems is 
the first one, why so? There are 3N — 6 internal degrees of freedom that span the 
global potential energy surface of an unconstrained iV-body system. Using for sim- 
plicity 10 discretization points per coordinate implies that of the order of 10 3N ~ 6 
electronic structure calculations are needed in order to map such a global potential 
energy surface. Thus, the computational workload for the first step grows roughly 
like ~ 10 N with increasing system size. This is what might be called the "dimen- 
sionality bottleneck" of calculations that rely on global potential energy surfaces, 
sec for instance the discussion on p. 420 in Ref. 254 . 

What is needed in ab initio molecular dynamics instead? Suppose that a useful 
trajectory consists of about 10 M molecular dynamics steps, i.e. 10 M electronic 
structure calculations arc needed to generate one trajectory. Furthermore, it is 
assumed that 10™ independent trajectories are necessary in order to average over 
different initial conditions so that \0 M + n ab initio molecular dynamics steps are 
required in total. Finally, it is assumed that each single point electronic structure 
calculation needed to devise the global potential energy surface and one ab initio 
molecular dynamics time step requires roughly the same amount of CPU time. Based 
on this truly simplistic order of magnitude estimate, the advantage of ab initio 
molecular dynamics vs. calculations relying on the computation of a global potential 
energy surface amounts to about io 3JV ~ 6 ~ M ~". The crucial point is that for a given 
statistical accuracy (that is for M and n fixed and independent on N) and for a 
given electronic structure method, the computational advantage of "on-the-fly" 
approaches grows like ^10^ with system size. 

Of course, considerable progress has been achieved in trajectory calculations by 
carefully selecting the discretization points and reducing their number, choosing so- 
phisticated representations and internal coordinates, exploiting symmetry etc. but 
basically the scaling ~ 10 N with the number of nuclei remains a problem. Other 
strategies consist for instance in reducing the number of active degrees of freedom 
by constraining certain internal coordinates, representing less important ones by a 
(harmonic) bath or friction, or building up the global potential energy surface in 
terms of few-body fragments. All these approaches, however, invoke approxima- 
tions beyond the ones of the electronic structure method itself. Finally, it is evident 
that the computational advantage of the "on-the-fly" approaches diminish as more 
and more trajectories are needed for a given (small) system. For instance extensive 
averaging over many different initial conditions is required in order to calculate 
quantitatively scattering or reactive cross sections. Summarizing this discussion, 
it can be concluded that ab initio molecular dynamics is the method of choice to 
investigate large and "chemically complex" systems. 

Quite a few review articles dealing with ab initio molecular dynamics appeared 

in the nineties 513,223,472,457,224,158,643,234,463,538,405 ^ ^ mterested reader ig fe _ 

ferred to them for various complementary viewpoints. In the present overview 
article, emphasis is put on both broadness of the approaches and depth of the pre- 
sentation. Concerning the broadness, the discussion starts from the Schrodingcr 
equation. Classical, Ehrenfest, Born-Oppcnhcimcr, and Car- Parrinello molecular 
dynamics are "derived" from the time-dependent mean-field approach that is ob- 
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tained after separating the nuclear and electronic degrees of freedom. The most 
extensive discussion is related to the features of the basic Car-Parrinello approach 
but all three ab initio approaches to molecular dynamics are contrasted and partly 
compared. The important issue of how to obtain the correct forces in these schemes 
is discussed in some depth. The most popular electronic structure theories imple- 
mented within ab initio molecular dynamics, density functional theory in the first 
place but also the Hartree-Fock approach, are sketched. Some attention is also 
given to another important ingredient in ab initio molecular dynamics, the choice 
of the basis set. 

Concerning the depth, the focus of the present discussion is clearly the im- 
plementation of both the basic Car-Parrincllo and Born-Oppenhcimcr molecular 
dynamics schemes in the CPMD package 142 . The electronic structure approach 
in CPMD is Hohenberg-Kohn-Sham density functional theory within a plane wave 
/ pseudopotential implementation and the Generalized Gradient Approximation. 
The formulae for energies, forces, stress, pseudopotentials, boundary conditions, 
optimization procedures, parallelization etc. are given for this particular choice to 
solve the electronic structure problem. One should, however, keep in mind that 
a variety of other powerful ab initio molecular dynamics codes are available (for 
instance CASTEP 116 , CP-PAW 143 , fhi98md 189 , NWChem 446 , VASP 663 ) which arc 
partly based on very similar techniques. The classic Car-Parrinello approach 108 
is then extended to other ensembles than the microcanonical one, other electronic 
states than the ground state, and to a fully quantum-mechanical representation of 
the nuclei. Finally, the wealth of problems that can be addressed using ab initio 
molecular dynamics is briefly sketched at the end, which also serves implicitly as 
the "Summary and Conclusions" section. 



2 Basic Techniques: Theory 

2.1 Deriving Classical Molecular Dynamics 

The starting point of the following discussion is non-relativistic quantum mechanics 
as formalized via the time-dependent Schrodinger equation 

i/iJU({rJ, {R 7 }; t) = W$({rJ, {R/}; t) (1) 
in its position representation in conjunction with the standard Hamiltonian 

& n2 V- & „2 V- e 2 c 2 Z/ x C 2 Z!Zj 



2Ml 1 ^2m c ^Z^\ Ti - Tj \ 2-,|R J _ r .| + 2 ; .|R J _R J | 

h 2 2 ^ h 2 

2Mj 1 V 2m e 

2Mj 



1~t 

^2 \l, ' ^ -In, . . . r , , , . , 

= ~ E^-E£-V 2 + K- c ({rJ, { R,}) 
i 

= -E^7 V ' + W e({r i },{R/}) (2) 

for the electronic {rj} and nuclear {R/} degrees of freedom. The more convenient 
atomic units (a.u.) will be introduced at a later stage for reasons that will soon 
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become clear. Thus, only the bare electron-electron, electron-nuclear, and nuclear- 
nuclear Coulomb interactions are taken into account. 

The goal of this section is to derive classical molecular dynamics 12 > 270 > 217 
starting from Schrodinger's wave equation and following the elegant route of 
Tully 650 > 651 . To this end, the nuclear and electronic contributions to the total 
wavefunction $({ri}, {R/}; t), which depends on both the nuclear and electronic 
coordinates, have to be separated. The simplest possible form is a product ansatz 

rt 

(3) 



$({r J, {R 7 }; t) « *({ ri }; t) x({R/}; t) exp 

L" J t 

where the nuclear and electronic wavefunctions are separately normalized to unity 
at every instant of time, i.e. (x\ t\x\ t) = 1 an d (^f;t\^f;t) = 1, respectively. In 
addition, a convenient phase factor 

E e = J drdK **({ ri }; t) X*({R/}; t) H e *({r J; t) X ({R,}; t) (4) 

was introduced at this stage such that the final equations will look nice; / drdR 
refers to the integration over all i = 1, . . . and I = 1, . . . variables {r^} and {R/}, 
respectively. It is mentioned in passing that this approximation is called a one- 
determinant or single-configuration ansatz for the total wavefunction, which at the 
end must lead to a mean-field description of the coupled dynamics. Note also that 
this product ansatz (excluding the phase factor) differs from the Born-Oppenheimer 
ansatz 340 > 350 for separating the fast and slow variables 

oo 

$Bo({ri}, {R/}; t) = *fe({rj, {R/})x/c({R/}; t) (5) 

k=0 

even in its one-determinant limit, where only a single electronic state k (evaluated 
for the nuclear configuration {R/}) is included in the expansion. 

Inserting the separation ansatz Eq. (3) into Eqs. (l)-(2) yields (after multiplying 
from the left by (^| and (x| and imposing energy conservation d(H) jdt = 0) the 
following relations 

^ = - E ^ v ?* + {/ dK x*« R '}; t)Vn-e({Ti}, {R/})x({R/}; *)} * (6) 

^ = -E^V?x+{/ dr**({r i };t)We({r i },{R/})*({ri};t)}x • (7) 

This set of coupled equations defines the basis of the time-dependent self-consistent 
field (TDSCF) method introduced as early as 1930 by Dirac 162 , see also Ref. 158 . 
Both electrons and nuclei move quantum-mechanically in time-dependent effective 
potentials (or self -consistently obtained average fields) obtained from appropriate 
averages (quantum mechanical expectation values (...)) over the other class of 
degrees of freedom (by using the nuclear and electronic wavefunctions, respectively) . 
Thus, the single-determinant ansatz Eq. (3) produces, as already anticipated, a 
mean field description of the coupled nuclear electronic quantum dynamics. This 
is the price to pay for the simplest possible separation of electronic and nuclear 
variables. 
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The next step in the derivation of classical molecular dynamics is the task to 
approximate the nuclei as classical point particles. How can this be achieved in the 
framework of the TDSCF approach, given one quantum-mechanical wave equa- 
tion describing all nuclei? A well -known route to extract classical mechanics from 
quantum mechanics in general starts with rewriting the corresponding wavefunction 

X({R / }; t) = A({Rj}; t) exp [iS({R/}; t)/h] (8) 

in terms of an amplitude factor A and a phase S which are both considered to be 
real and A)0 in this polar representation, see for instance Refs. 16 3,425,535 After 
transforming the nuclear wavefunction in Eq. (7) accordingly and after separating 
the real and imaginary parts, the TDSCF equation for the nuclei 



2 A 
A 



(9) 



f + E^(v^) 2 +/^ C *^£^I 

i J i 

^^MM + i;^W|.o do) 

is (exactly) re-expressed in terms of the new variables A and S. This so-called 
"quantum fluid dynamical representation" Eqs. (9)-(10) can actually be used to 
solve the time-dependent Schrodingcr equation 160 . The relation for A, Eq. (10), 
can be rewritten as a continuity equation 16 M25,535 the h e }p G f the identi- 

fication of the nuclear density \x\ 2 = A 2 as directly obtained from the definition 
Eq. (8). This continuity equation is independent of h and ensures locally the con- 
servation of the particle probability \x\ 2 associated to the nuclei in the presence of 
a flux. 

More important for the present purpose is a more detailed discussion of the 
relation for S, Eq. (9). This equation contains one term that depends on h, a 
contribution that vanishes if the classical limit 

f + E^(^) 2 + /*W = o (ID 

is taken as H — > 0; an expansion in terms of H would lead to a hierarchy of semi- 
classical methods 425 - 259 . The resulting equation is now isomorphic to equations of 
motion in the Hamilton- J acobi formulation 244 ' 540 

— +W({R / },{V / 5})=0 (12) 

of classical mechanics with the classical Hamilton function 

W({R,}, {P,}) = T({Pj}) + y({R 7 }) (13) 

defined in terms of (generalized) coordinates {R/} and their conjugate momenta 
{P/}. With the help of the connecting transformation 

Pi = V/S (14) 



7 



the Newtonian equation of motion Pj = — V/V({Rj}) corresponding to Eq. (11) 

dPi 



= -Vj / dr **W C * or 

M 7 R 7 (t) = -V/ y drtf*W e tf (15) 

= -V/K E ({R/(t)}) (16) 

can be read off. Thus, the nuclei move according to classical mechanics in an 
effective potential due to the electrons. This potential is a function of only the 
nuclear positions at time t as a result of averaging H c over the electronic degrees 
of freedom, i.e. computing its quantum expectation value (WlWcl^), while keeping 
the nuclear positions fixed at their instantaneous values {R/(t)}. 

However, the nuclear wavefunction still occurs in the TDSCF equation for the 
electronic degrees of freedom and has to be replaced by the positions of the nuclei for 
consistency. In this case the classical reduction can be achieved simply by replacing 
the nuclear density |x({R/}; t)\ 2 in Eq. (6) in the limit h — > by a product of delta 
functions Yli ^(Rj — R/(^)) centered at the instantaneous positions {R/(i)} of the 
classical nuclei as given by Eq. (15). This yields e.g. for the position operator 

ydR x *({R/};*)R/x({R/};*) R/(t) (17) 

the required expectation value. This classical limit leads to a time dependent wave 
equation for the electrons 

ih ~dT = ~ E ^ V ?* + K-e({rJ, {R,(t)»* 

i 

= H e ({r i },{K I (t)}) *({rJ,{R 7 };i) (18) 

which evolve self consistently as the classical nuclei are propagated via Eq. (15). 
Note that now H e and thus \t depend parametrically on the classical nuclear posi- 
tions {R/(t)} at time t through Vn_ e ({r,}, {Rj(t)}). This means that feedback 
between the classical and quantum degrees of freedom is incorporated in both 
directions (at variance with the "classical path" or Mott non-SCF approach to 
dynamics 650 < 651 ). 

The approach relying on solving Eq. (15) together with Eq. (18) is sometimes 
called "Ehrenfest molecular dynamics" in honor of Ehrcnfest who was the first to 
address the question a of how Newtonian classical dynamics can be derived from 
Schrodinger's wave equation 174 . In the present case this leads to a hybrid or 
mixed approach because only the nuclei are forced to behave like classical particles, 
whereas the electrons are still treated as quantum objects. 

Although the TDSCF approach underlying Ehrcnfest molecular dynamics 
clearly is a mean-field theory, transitions between electronic states are included 



"The opening statement of Ehrcnfest's famous 1927 paper 174 reads: 

"Es ist wiinschenswert, die folgcnde Frage moglichst clemcntar bcantwortcn zu konncn: Welcher 
Riickblick ergibt sich vom Standpunkt der Quantenmechanik auf die Newtonschen Grundgleichun- 
gen der klassischen Mechanik?" 
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in this scheme. This can be made evident by expanding the electronic wavcfunc- 
tion ^/ (as opposed to the total wavefunction $ according to Eq. (5)) in terms of 
many electronic states or determinants 

oo 

f({r i}, {R/}; t) = J2 Cfc(t)*fc({ri}; {R/}) (19) 

fe=0 

with complex coefficients {cfc(t)}. In this case, the coefficients {|cfc(i)| 2 } (with 
Sfc l c fc(^)| 2 = 1) describe explicitly the time evolution of the populations (occupa- 
tions) of the different states {k} whereas interferences are included via the {c^ci^k} 
contributions. One possible choice for the basis functions {^fe} is the adiabatic basis 
obtained from solving the time -independent electronic Schrodinger equation 

We({r j };{Rj})*fc = Sfc({Rj})*fc({r i };{Rj}) , (20) 

where {R/} are the instantaneous nuclear positions at time t according to Eq. (15). 
The actual equations of motion in terms of the expansion coefficients {ck} are 
presented in Sect. 2.2. 

At this stage a further simplification can be invoked by restricting the total 
electronic wave function ^ to be the ground state wave function of Tt e at each 
instant of time according to Eq. (20) and |co(i)| 2 = 1 in Eq. (19). This should be a 
good approximation if the energy difference between *S? and the first excited state 

is everywhere large compared to the thermal energy k^T , roughly speaking. In 
this limit the nuclei move according to Eq. (15) on a single potential energy surface 

Vf = Jdr = £o({R/}) (21) 

that can be computed by solving the time-independent electronic Schrodinger equa- 
tion Eq. (20) 

H c y = E y , (22) 

for the ground state only. This leads to the identification V® = E via Eq. (21), 
i.e. in this limit the Ehrenfest potential is identical to the ground-state Born 
Oppcnhcimcr potential. 

As a consequence of this observation, it is conceivable to decouple the task of 
generating the nuclear dynamics from the task of computing the potential energy 
surface. In a first step Eq is computed for many nuclear configurations by solving 
Eq. (22). In a second step, these data points are fitted to an analytical functional 
form to yield a global potential energy surface 539 , from which the gradients can be 
obtained analytically. In a third step, the Newtonian equation of motion Eq. (16) 
is solved on this surface for many different initial conditions, producing a "swarm" 
of classical trajectories. This is, in a nutshell, the basis of classical trajectory cal- 
culations on global potential energy surfaces 649 ' 541 . 

As already alluded to in the general introduction, such approaches suffer severely 
from the "dimensionality bottleneck" as the number of active nuclear degrees of 
freedom increases. One traditional way out of this dilemma is to approximate the 
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global potential energy surface 



yE „ K approx ({R/}) = £ „ l(Rj) + £ ^ (Rj ; Rj) 



1=1 I(J 



N 



+ «3(Rj,Rj,R«-) + -- 



(23) 



in terms of a truncated expansion of many-body contributions 253 > 12 > 270 . At this 
stage, the electronic degrees of freedom are replaced by interaction potentials {v n } 
and are not featured as explicit degrees of freedom in the equations of motion. Thus, 
the mixed quantum / classical problem is reduced to purely classical mechanics, 
once the {v n } are determined. Classical molecular dynamics 



relies crucially on this idea, where typically only two-body V2 or three-body V3 
interactions are taken into account 12 > 270 ; although more sophisticated models to 
include non- additive interactions such as polarization exist. This amounts to a 
dramatic simplification and removes the dimensionality bottleneck as the global 
potential surface is constructed from a manageable sum of additive few-body con- 
tributions — at the price of introducing a drastic approximation and of basically 
excluding chemical transformations from the realm of simulations. 

As a result of this derivation, the essential assumptions underlying classical 
molecular dynamics become transparent: the electrons follow adiabatically the clas- 
sical nuclear motion and can be integrated out so that the nuclei evolve on a single 
Born-Oppenheimer potential energy surface (typically but not necessarily given by 
the electronic ground state), which is in general approximated in terms of few-body 
interactions. 

Actually, classical molecular dynamics for many-body systems is only made 
possible by somehow decomposing the global potential energy. In order to illustrate 
this point consider the simulation of N = 500 Argon atoms in the liquid phase 175 
where the interactions can faithfully be described by additive two-body terms, i.e. 
^approx^p^}.) _ j2f {J v 2 (|Rj-Rj|). Thus, the determination of the pair potential 
v 2 from ah initio electronic structure calculations amounts to computing and fitting 
a one dimensional function. The corresponding task to determine a global potential 
energy surface amounts to doing that in about 10 1500 dimensions, which is simply 
impossible (and on top of that not necessary for Nobel gases!). 

2.2 Ehrenfest Molecular Dynamics 

A way out of the dimensionality bottleneck other than to approximate the global 
potential energy surface Eq. (23) or to reduce the number of active degrees of free- 
dom is to take seriously the classical nuclei approximation to the TDSCF equations, 
Eq. (15) and (18). This amounts to computing the Ehrenfest force by actually solv- 



MjB,(t) 



V/V e approx ({Ri(t)}) 



(24) 
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ing numerically 



MjR/(t) 



= -V/ J drV 



-V/(*|W e |*> 

-V/ (H c ) 
7 /l 



(25) 



K 



^2^T Vl2 + K - c({rJ ' {R/(<)}) 



(26) 



the coupled set of equations simultaneously. Thereby the a priori construction 
of any type of potential energy surface is avoided from the outset by solving the 
time-dependent electronic Schrodinger equation "on-the-fly" . This allows one to 
compute the force from Vi(H e ) for each configuration {Rj(t)} generated by molec- 
ular dynamics; see Sect. 2.5 for the issue of using the so-called "Hcllmann Fcynman 
forces" instead. 

The corresponding equations of motion in terms of the adiabatic basis Eq. (20) 
and the time-dependent expansion coefficients Eq. (19) read 650 > 651 



M&iit) = -J2 \ck{t)\ 2 V!E k - 4d (E k - Ei) d 
k k.i 

ihc k (t) = c k (t)E k - ihy^ci(t)Kid 



if 



where the coupling terms are given by 

df({Rj(t)}) = J dr*J;Vi*« 



(27) 
(28) 

(29) 



with the property dj fc = 0. The Ehrenfest approach is thus seen to include rigor- 
ously non-adiabatic transitions between different electronic states ^ k and within 
the framework of classical nuclear motion and the mean-field (TDSCF) approxi- 
mation to the electronic structure, see e.g. Refs. 650 > 651 for reviews and for instance 
Ref. 532 for an implementation in terms of time-dependent density functional the- 
ory. 

The restriction to one electronic state in the expansion Eq. (19), which is in 
most cases the ground state 'to, leads to 



MjRi(i) = -V/ (#„ |«e|*0> 



th- 



at 



(30) 
(31) 



as a special case of Eqs. (25)-(26); note that Tt e is time-dependent via the nuclear 
coordinates {R/(t)}. A point worth mentioning here is that the propagation of the 
wavefunction is unitary, i.e. the wavefunction preserves its norm and the set of 
orbitals used to build up the wavefunction will stay orthonormal, see Sect. 2.6. 
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Ehrenfest molecular dynamics is certainly the oldest approach to "on-the-fly" 
molecular dynamics and is typically used for collision- and scattering-type prob- 
lems 154 i 649 > 426 ! 532 . However, it was never in widespread use for systems with many 
active degrees of freedom typical for condensed matter problems for reasons that 
will be outlined in Sec. 2.6 (although a few exceptions exist 553 < 34 > 203 > 617 but here 
the number of explicitly treated electrons is fairly limited with the exception of 
Rcf. 617 ). 

2.3 Born-Oppenheimer Molecular Dynamics 

An alternative approach to include the electronic structure in molecular dynamics 
simulations consists in straightforwardly solving the static electronic structure prob- 
lem in each molecular dynamics step given the set of fixed nuclear positions at that 
instance of time. Thus, the electronic structure part is reduced to solving a time- 
independent quantum problem, e.g. by solving the time-independent Schrodingcr 
equation, concurrently to propagating the nuclei via classical molecular dynamics. 
Thus, the time-dependence of the electronic structure is a consequence of nuclear 
motion, and not intrinsic as in Ehrenfest molecular dynamics. The resulting Born 
Oppcnhcimcr molecular dynamics method is defined by 

M/Rj(t) = -Vjmin{<tf |«e|*o>} (32) 

Eo"£ = H c *o (33) 

for the electronic ground state. A deep difference with respect to Ehrenfest dy- 
namics concerning the nuclear equation of motion is that the minimum of (Tie) 
has to be reached in each Born-Oppenheimer molecular dynamics step according 
to Eq. (32). In Ehrenfest dynamics, on the other hand, a wavefunction that min- 
imized (H c ) initially will also stay in its respective minimum as the nuclei move 
according to Eq. (30)! 

A natural and straightforward extension 281 of ground-state Born-Oppenheimer 
dynamics is to apply the same scheme to any excited electronic state without 
considering any interferences. In particular, this means that also the "diagonal 
correction terms" 340 

r>* fc ({Rj(t)}) = ~Jdr *£V?* fc (34) 

arc always neglected; the inclusion of such terms is discussed for instance in 
Refs. 650 > 651 . These terms renormalize the Born-Oppenheimer or "clamped nu- 
clei" potential energy surface Ek of a given state \Pfc (which might also be the 
ground state \fo) and lead to the so-called "adiabatic potential energy surface" 
of that state 340 . Whence, Born-Oppenheimer molecular dynamics should not be 
called "adiabatic molecular dynamics" , as is sometime done. 

It is useful for the sake of later reference to formulate the Born-Oppenheimer 
equations of motion for the special case of effective one-particle Hamiltonians. This 
might be the Hartree-Fock approximation defined to be the variational minimum 
of the energy expectation value \He \ *o) given a single Slater determinant "Jo = 
det-j^i} subject to the constraint that the one-particle orbitals ipi arc orthonormal 
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{ipi \tpj ) — Sij. The corresponding constraint minimization of the total energy with 
respect to the orbitals 



min{(* | Wc |v]/ )} 

{>/>>} 



(35) 

{<Vi|^>=5«} 



can be cast into Lagrange's formalism 

c = - <* \n c \ * ) + A ^ 1^- ) - 5 «) ( 36 ) 

«>.J 

where are the associated Lagrangian multipliers. Unconstrained variation of 
this Lagrangian with respect to the orbitals 

leads to the well-known Hartree-Fock equations 

w e HF v>i = E A ^ ( 38 ) 

3 

as derived in standard text books 604 > 418 ; the diagonal canonical form H^ipi = t^i 
is obtained after a unitary transformation and H^ F denotes the effective one 
particle Hamiltonian, see Sect. 2.7 for more details. The equations of motion 
corresponding to Eqs. (32)-(33) read 

Mj-Riit) = -V/min{(* |^ F |*o)} (39) 
{</>;} 

= -^+^A^ (40) 

3 

for the Hartree-Fock case. A similar set of equations is obtained if Hohcnberg- 
Kohn-Sham density functional theory 458 ' 168 is used, where Tt^} F has to be replaced 
by the Kohn-Sham effective one-particle Hamiltonian H^ s , see Sect. 2.7 for more 
details. Instead of diagonalizing the one-particle Hamiltonian an alternative but 
equivalent approach consists in directly performing the constraint minimization 
according to Eq. (35) via nonlinear optimization techniques. 

Early applications of Born-Oppcnhcimcr molecular dynamics were performed 
in the framework of a scmicmpirical approximation to the electronic structure prob- 
lem 669 ' 671 . But only a few years later an ab initio approach was implemented within 
the Hartree-Fock approximation 365 . Born-Oppcnhcimcr dynamics started to be- 
come popular in the early nineties with the availability of more efficient electronic 
structure codes in conjunction with sufficient computer power to solve "interesting 
problems", see for instance the compilation of such studies in Table 1 in a recent 
overview article 82 . 

Undoubtedly, the breakthrough of Hohenberg-Kohn-Sham density functional 
theory in the realm of chemistry - which took place around the same time - also 
helped a lot by greatly improving the "price / performance ratio" of the electronic 
structure part, see e.g. Refs. 694 > 590 . A third and possibly the crucial reason that 
boosted the field of ab initio molecular dynamics was the pioneering introduction of 
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the Car-Parrinello approach 108 , see also Fig. 1. This technique opened novel av- 
enues to treat large-scale problems via ah initio molecular dynamics and catalyzed 
the entire field by making "interesting calculations" possible, see also the closing 
section on applications. 

2.4 Car-Parrinello Molecular Dynamics 
2.4-1 Motivation 

A non-obvious approach to cut down the computational expenses of molecular dy- 
namics which includes the electrons in a single state was proposed by Car and 
Parrinello in 1985 108 . In retrospect it can be considered to combine the advan- 
tages of both Ehrenfest and Born-Oppcnhcimcr molecular dynamics. In Ehrcnfcst 
dynamics the time scale and thus the time step to integrate Eqs. (30) and (31) 
simultaneously is dictated by the intrinsic dynamics of the electrons. Since elec- 
tronic motion is much faster than nuclear motion, the largest possible time step 
is that which allows to integrate the electronic equations of motion. Contrary 
to that, there is no electron dynamics whatsoever involved in solving the Born- 
Oppenheimer Eqs. (32)-(33), i.e. they can be integrated on the time scale given 
by nuclear motion. However, this means that the electronic structure problem 
has to be solved self-consistently at each molecular dynamics step, whereas this is 
avoided in Ehrenfest dynamics due to the possibility to propagate the wavefunc- 
tion by applying the Hamiltonian to an initial wavefunction (obtained e.g. by one 
self-consistent diagonalization) . 

From an algorithmic point of view the main task achieved in ground-state 
Ehrenfest dynamics is simply to keep the wavefunction automatically minimized 
as the nuclei are propagated. This, however, might be achieved - in principle - by 
another sort of deterministic dynamics than first-order Schrodingcr dynamics. In 
summary, the "Best of all Worlds Method" should (i) integrate the equations of 
motion on the (long) time scale set by the nuclear motion but nevertheless (ii) take 
intrinsically advantage of the smooth time evolution of the dynamically evolving 
electronic subsystem as much as possible. The second point allows to circumvent 
explicit diagonalization or minimization to solve the electronic structure problem 
for the next molecular dynamics step. Car-Parrinello molecular dynamics is an ef- 
ficient method to satisfy requirement (ii) in a numerically stable fashion and makes 
an acceptable compromise concerning the length of the time step (i). 

2-4-2 Car-Parrinello Lagrangian and Equations of Motion 

The basic idea of the Car-Parrinello approach can be viewed to exploit the 
quantum-mechanical adiabatic time-scale separation of fast electronic and slow 
nuclear motion by transforming that into classical-mechanical adiabatic energy- 
scale separation in the framework of dynamical systems theory. In order to achieve 
this goal the two-component quantum / classical problem is mapped onto a two- 
component purely classical problem with two separate energy scales at the expense 
of loosing the explicit time-dependence of the quantum subsystem dynamics. Fur- 
thermore, the central quantity, the energy of the electronic subsystem (^ol^ol^o) 
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evaluated with some wavefunction ^o, is certainly a function of the nuclear posi- 
tions {R/}. But at the same time it can be considered to be a functional of the 
wavefunction Wo and thus of a set of one particle orbitals {ipi} (or in general of 
other functions such as two-particle gcminals) used to build up this wavefunction 
(being for instance a Slater determinant Wo = det{i/>j} or a combination thereof). 
Now, in classical mechanics the force on the nuclei is obtained from the deriva- 
tive of a Lagrangian with respect to the nuclear positions. This suggests that a 
functional derivative with respect to the orbitals, which are interpreted as classical 
fields, might yield the force on the orbitals, given a suitable Lagrangian. In addi- 
tion, possible constraints within the set of orbitals have to be imposed, such as e.g. 
orthonormality (or generalized orthonormality conditions that include an overlap 
matrix) . 

Car and Parrincllo postulated the following class of Lagrangians 108 

£ CP = ^^M/Rf + ^^^ (^olWcl^o) + constraints, (41) 

^ j ' potential energy orthonormality 

kinetic energy 

to serve this purpose. The corresponding Newtonian equations of motion are ob- 
tained from the associated Euler-Lagrange equations 



dt 9R / <9R/ 
±^L = ^L (43 ) 

like in classical mechanics, but here for both the nuclear positions and the orbitals; 
note = (ipi\ and that the constraints are holonomic 244 . Following this route of 
ideas, generic Car-Parrinello equations of motion are found to be of the form 

M 7 R 7 (t) = -— - (* |H e |*o) + 75=- {constraints} (44) 

s s 

Mi^M*) = -— (*o|^o|*o) + T7- {constraints} (45) 

where \ii (= fi) are the "fictitious masses" or inertia parameters assigned to the 
orbital degrees of freedom; the units of the mass parameter [i are energy times a 
squared time for reasons of dimensionality. Note that the constraints within the 
total wavefunction lead to "constraint forces" in the equations of motion. Note also 
that these constraints 

constraints — constraints ({^},{R/}) (46) 

might be a function of both the set of orbitals {ipi} and the nuclear positions {R/}. 
These dependencies have to be taken into account properly in deriving the Car- 
Parrinello equations following from Eq. (41) using Eqs. (42)-(43), see Sect. 2.5 for 
a general discussion and see e.g. Ref. 351 for a case with an additional dependence 
of the wavefunction constraint on nuclear positions. 
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According to the Car-Parrinello equations of motion, the nuclei evolve in time 
at a certain (instantaneous) physical temperature ex MrR-j, whereas a "fic- 
titious temperature" oc [ii(ipi\ipi) is associated to the electronic degrees of 
freedom. In this terminology, "low electronic temperature" or "cold electrons" 
means that the electronic subsystem is close to its instantaneous minimum energy 
min{^.}(\E'o|^c| v I'o), i.e. close to the exact Born-Oppenheimer surface. Thus, a 
ground-state wavefunction optimized for the initial configuration of the nuclei will 
stay close to its ground state also during time evolution if it is kept at a sufficiently 
low temperature. 

The remaining task is to separate in practice nuclear and electronic motion such 
that the fast electronic subsystem stays cold also for long times but still follows 
the slow nuclear motion adiabatically (or instantaneously). Simultaneously, the 
nuclei are nevertheless kept at a much higher temperature. This can be achieved 
in nonlinear classical dynamics via decoupling of the two subsystems and (quasi-) 
adiabatic time evolution. This is possible if the power spectra stemming from 
both dynamics do not have substantial overlap in the frequency domain so that 
energy transfer from the "hot nuclei" to the "cold electrons" becomes practically 
impossible on the relevant time scales. This amounts in other words to imposing and 
maintaining a metastability condition in a complex dynamical system for sufficiently 
long times. How and to which extend this is possible in practice was investigated in 
detail in an important investigation based on well-controlled model systems 467 > 468 
(see also Sects. 3.2 and 3.3 in Ref. 513 ), with more mathematical rigor in Rcf. 86 , 
and in terms of a generalization to a second level of adiabaticity in Ref. 411 . 



2.4.3 Why Does the Car-Parrinello Method Work ? 

In order to shed light on the title question, the dynamics generated by the Car- 
Parrinello Lagrangian Eq. (41) is analyzed in more detail invoking a "classical 
dynamics perspective" of a simple model system (eight silicon atoms forming a 
periodic diamond lattice, local density approximation to density functional theory, 
normconserving pseudopotentials for core electrons, plane wave basis for valence 
orbitals, 0.3 fs time step with \i = 300 a.u., in total 20 000 time steps or 6.3 ps, 
for full details see Ref. 467 ); a concise presentation of similar ideas can be found 
in Ref. 110 . For this system the vibrational density of states or power spectrum 
of the electronic degrees of freedom, i.e. the Fourier transform of the statistically 
averaged velocity autocorrelation function of the classical fields 



/(<*>) = f Q dtcos(Lut)J2(^t\^0) (47) 



is compared to the highest-frequency phonon mode w™ ax of the nuclear subsystem 
in Fig. 2. From this figure it is evident that for the chosen parameters the nuclear 
and electronic subsystems are dynamically separated: their power spectra do not 
overlap so that energy transfer from the hot to the cold subsystem is expected to 
be prohibitively slow, see Sect. 3.3 in Ref. 513 for a similar argument. 

This is indeed the case as can be verified in Fig. 3 where the conserved energy 
S C ons, physical total energy -E p hy S > electronic energy V c , and fictitious kinetic energy 
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Figure 2. Vibrational density of states Eq. (47) (continuous spectrum in upper part) and harmonic 
approximation thereof Eq. (52) (stick spectrum in lower part) of the electronic degrees of freedom 
compared to the highest-frequency phonon mode o;™ ax (triangle) for a model system; for further 
details see text. Adapted from Ref. 467 . 



of the electrons T c 



E cons = \»i & ) + _C l M ^ 2 I + (*o|«e|*0> 



£phys = \ M ^ 2 I + (*o|«e|*0> = Scons - T c 



V c = (*o|H e |*o) 

T c = (jPi 



(48) 

(49) 

(50) 
(51) 



are shown for the same system as a function of time. First of all, there should be a 
conserved energy quantity according to classical dynamics since the constraints are 
holonomic 244 . Indeed "the Hamiltonian" or conserved energy E cons is a constant of 
motion (with relative variations smaller than 10 -6 and with no drift), which serves 
as an extremely sensitive check of the molecular dynamics algorithm. Contrary 
to that the electronic energy V c displays a simple oscillation pattern due to the 
simplicity of the phonon modes. 

Most importantly, the fictitious kinetic energy of the electrons T e is found to 
perform bound oscillations around a constant, i.e. the electrons "do not heat up" 
systematically in the presence of the hot nuclei; note that T c is a measure for devi- 
ations from the exact Born-Oppenhcimcr surface. Closer inspection shows actually 
two time scales of oscillations: the one visible in Fig. 3 stems from the drag exerted 
by the moving nuclei on the electrons and is the mirror image of the V e fluctuations. 
Superimposed on top of that (not shown, but see Fig. 4(b)) arc small-amplitude 
high frequency oscillations intrinsic to the fictitious electron dynamics with a period 
of only a fraction of the visible mode. These oscillations are actually instrumental 
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Figure 3. Various energies Eqs. (48)— (51) for a model system propagated via Car-Parrinello molec- 
ular dynamics for at short (up to 300 fs), intermediate, and long times (up to 6.3 ps); for further 
details see text. Adapted from Ref. 467 . 



for the stability of the Car-Parrinello dynamics, vide infra. But already the visible 
variations are three orders of magnitude smaller than the physically meaningful os- 
cillations of V c . As a result, -E p h ys defined as E cons — T c or cquivalcntly as the sum 
of the nuclear kinetic energy and the electronic total energy (which serves as the 
potential energy for the nuclei) is essentially constant on the relevant energy and 
time scales. Thus, it behaves approximately like the strictly conserved total energy 
in classical molecular dynamics (with only nuclei as dynamical degrees of freedom) 
or in Born Oppcnhcimcr molecular dynamics (with fully optimized electronic de- 
grees of freedom) and is therefore often denoted as the "physical total energy" . 
This implies that the resulting physically significant dynamics of the nuclei yields 
an excellent approximation to microcanonical dynamics (and assuming ergodicity 
to the microcanonical ensemble). Note that a different explanation was advocated 
in Ref. 470 (see also Ref. 472 , in particular Sect. VIII. B and C), which was however 
revised in Ref. no . A discussion similar in spirit to the one outlined here 467 is 
provided in Ref. 513 , see in particular Sect. 3.2 and 3.3. 

Given the adiabatic separation and the stability of the propagation, the central 
question remains if the forces acting on the nuclei are actually the "correct" ones 
in Car-Parrinello molecular dynamics. As a reference serve the forces obtained 
from full self-consistent minimizations of the electronic energy min^.} (^ol^ol^o) 
at each time step, i.e. Born-Oppenheimer molecular dynamics with extremely well 
converged wavefunctions. This is indeed the case as demonstrated in Fig. 4(a): 
the physically meaningful dynamics of the ^-component of the force acting on one 
silicon atom in the model system obtained from stable Car-Parrinello fictitious 
dynamics propagation of the electrons and from iterative minimizations of the elec- 
tronic energy are extremely close. 

Better resolution of one oscillation period in (b) reveals that the gross devia- 
tions are also oscillatory but that they arc four orders of magnitudes smaller than 
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Figure 4. (a) Comparison of the x— component of the force acting on one atom of a model system 
obtained from Car-Parrinello (solid line) and well-converged Born-Oppcnhcimcr (dots) molecular 
dynamics, (b) Enlarged view of the difference between Car-Parrinello and Born-Oppenheimer 
forces; for further details see text. Adapted from Ref. 467 . 



the physical variations of the force resolved in Fig. 4(a). These correspond to the 
"large amplitude" oscillations of T c visible in Fig. 3 due to the drag of the nuclei 
exerted on the quasi adiabatically following electrons having a finite dynamical 
mass /z. Note that the inertia of the electrons also dampens artificially the nuclear 
motion (typically on a few percent scale, see Sect. V.C.2 in Ref. 75 for an anal- 
ysis and a renormalization correction of Mj) but decreases as the fictitious mass 
approaches the adiabatic limit \i — > 0. Superimposed on the gross variation in (b) 
are again high-frequency bound oscillatory small-amplitude fluctuations like for T c . 
They lead on physically relevant time scales (i.e. those visible in Fig. 4(a)) to "av- 
eraged forces" that are very close to the exact ground-state Born-Oppenheimer 
forces. This feature is an important ingredient in the derivation of adiabatic dy- 
namics 467 < 411 . 

In conclusion, the Car-Parrinello force can be said to deviate at most instants of 
time from the exact Born-Oppenheimer force. However, this does not disturb the 
physical time evolution due to (i) the smallness and boundedness of this difference 
and (ii) the intrinsic averaging effect of small-amplitude high-frequency oscillations 
within a few molecular dynamics time steps, i.e. on the sub-femtosecond time scale 
which is irrelevant for nuclear dynamics. 



2.4-4 How to Control Adiabaticity ? 

An important question is under which circumstances the adiabatic separation can 
be achieved, and how it can be controlled. A simple harmonic analysis of the 
frequency spectrum of the orbital classical fields close to the minimum defining the 
ground state yields 467 
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where €j and e, are the eigenvalues of occupied and unoccupied orbitals, respec- 
tively; see Eq. (26) in Rcf. 467 for the case where both orbitals are occupied ones. 
It can be seen from Fig. 2 that the harmonic approximation works faithfully as 
compared to the exact spectrum; see Ref. 471 and Sect. IV. A in Ref. 472 for a more 
general analysis of the associated equations of motion. Since this is in particu- 
lar true for the lowest frequency w™ m , the handy analytic estimate for the lowest 
possible electronic frequency 



shows that this frequency increases like the square root of the electronic energy 
difference E gap between the lowest unoccupied and the highest occupied orbital. 
On the other hand it increases similarly for a decreasing fictitious mass parameter 
/i. 

In order to guarantee the adiabatic separation, the frequency difference w™ m — 
should be large, see Sect. 3.3 in Ref. 513 for a similar argument. But both 
the highest phonon frequency w™ ax and the energy gap E gap are quantities that a 
dictated by the physics of the system. Whence, the only parameter in our hands 
to control adiabatic separation is the fictitious mass, which is therefore also called 
"adiabaticity parameter". However, decreasing ix not only shifts the electronic 
spectrum upwards on the frequency scale, but also stretches the entire frequency 
spectrum according to Eq. (52) . This leads to an increase of the maximum frequency 
according to 



where _E cut is the largest kinetic energy in an expansion of the wavefunction in 
terms of a plane wave basis set, see Sect. 3.1.3. 

At this place a limitation to decrease /x arbitrarily kicks in due to the maximum 
length of the molecular dynamics time step At max that can be used. The time step 
is inversely proportional to the highest frequency in the system, which is o>™ ax and 
thus the relation 



governs the largest time step that is possible. As a consequence, Car-Parrinello 
simulators have to find their way between Scylla and Charybdis and have to make 
a compromise on the control parameter /x; typical values for large-gap systems are 
/x = 500-1500 a.u. together with a time step of about 5-10 a.u. (0.12-0.24 fs). 
Recently, an algorithm was devised that optimizes /x during a particular simulation 
given a fixed accuracy criterion 87 . Note that a poor man's way to keep the time 
step large and still increase \i in order to satisfy adiabaticity is to choose heavier 
nuclear masses. That depresses the largest phonon or vibrational frequency w™ ax 
of the nuclei (at the cost of renormalizing all dynamical quantities in the sense of 
classical isotope effects). 




(53) 




(54) 




(55) 
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Up to this point the entire discussion of the stability and adiabaticity issues 
was based on model systems, approximate and mostly qualitative in nature. But 
recently it was actually proven 86 that the deviation or the absolute error A M of the 
Car-Parrinello trajectory relative to the trajectory obtained on the exact Born 
Oppenheimer potential energy surface is controlled by ^t: 
Theorem 1 iv.): There are constants C)0 and [i*)0 such that 



An = |B/*(t) -R°(t)| + || V";*) ~ I^V)! < Cy /2 , 0<t<T (56) 



for all values of the parameter \x satisfying 0(/i < [/,*, where up to time T)0 there 
exists a unique nuclear trajectory on the exact Born-Oppenhcimcr surface with 
w™ ln )0 for < t < T, i.e. there is "always" a finite electronic excitation gap. 
Here, the superscript /i or indicates that the trajectory was obtained via Car- 
Parrinello molecular dynamics using a finite mass fx or via dynamics on the exact 
Born-Oppenheimer surface, respectively. Note that not only the nuclear trajectory 
is shown to be close to the correct one, but also the wavefunction is proven to stay 
close to the fully converged one up to time T. Furthermore, it was also investi- 
gated what happens if the initial wavefunction at t = is not the minimum of the 
electronic energy (H c ) but trapped in an excited state. In this case it is found that 
the propagated wavefunction will keep on oscillating at t}0 also for ji — > and not 
even time averages converge to any of the cigenstates. Note that this does not pre- 
clude Car-Parrinello molecular dynamics in excited states, which is possible given 
a properly "minimizable" expression for the electronic energy, see e.g. Refs. 281 > 214 . 
However, this finding might have crucial implications for electronic level-crossing 
situations. 

What happens if the electronic gap is very small or even vanishes E gap — > 
as is the case for metallic systems? In this limit, all the above-given arguments 
break down due to the occurrence of zero frequency electronic modes in the power 
spectrum according to Eq. (53), which necessarily overlap with the phonon spec- 
trum. Following an idea of Sprik 583 applied in a classical context it was shown 
that the coupling of separate Nose-Hoover thermostats 12 > 270 > 217 to the nuclear and 
electronic subsystem can maintain adiabaticity by counterbalancing the energy flow 
from ions to electrons so that the electrons stay "cool" 74 ; see Rcf. 204 for a simi- 
lar idea to restore adiabaticity. Although this method is demonstrated to work in 
practice , this ad hoc cure is not entirely satisfactory from both a theoretical and 
practical point of view so that the well - controlled Born-Oppenheimer approach is 
recommended for strongly metallic systems. An additional advantage for metal- 
lic systems is that the latter is also better suited to sample many k-points (see 
Sect. 3.1.3), allows easily for fractional occupation numbers 458 ' 168 ; and can handle 
efficiently the so-called charge sloshing problem 472 . 



and the fictitious kinetic energy satisfies 




< t < T 



(57) 
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2.4-5 The Quantum Chemistry Viewpoint 



In order to understand Car-Parrinello molecular dynamics also from the "quantum 
chemistry perspective" , it is useful to formulate it for the special case of the Hartree- 
Fock approximation using 



£ CP = E \ m i&i + E \^ 



<*o|H c HF |*o) + £ ((^ \^)~ Sij) ■ (58) 

hi 

The resulting equations of motion 

AfjRi(t) = -V/(*o|^ F |*o) (59) 

IH$i{t) = -H? F ^ + E AyVj (60) 
j 

are very close to those obtained for Born -Oppcnhcimcr molecular dynamics 
Eqs. (39)-(40) except for (i) no need to minimize the electronic total energy ex- 
pression and (ii) featuring the additional fictitious kinetic energy term associated 
to the orbital degrees of freedom. It is suggestive to argue that both sets of equa- 
tions become identical if the term \ni"<pi{t)\ is small at any time t compared to the 
physically relevant forces on the right-hand-side of both Eq. (59) and Eq. (60). 
This term being zero (or small) means that one is at (or close to) the minimum of 
the electronic energy |7Y^ IF | \I/o) since time derivatives of the orbitals {tpi} can 
be considered as variations of Wo an d thus of the expectation value (7i,P F ) itself. 
In other words, no forces act on the wavefunction if Hiipi = 0. In conclusion, the 
Car-Parrinello equations are expected to produce the correct dynamics and thus 
physical trajectories in the microcanonical ensemble in this idealized limit. But 
if \/iiipi(t)\ is small for all i, this also implies that the associated kinetic energy 

= J2i A t i(V'i|V , i)/2 is small, which connects these more qualitative arguments 
with the previous discussion 467 . 

At this stage, it is also interesting to compare the structure of the Lagrangian 
Eq. (58) and the Euler-Lagrange equation Eq. (43) for Car-Parrinello dynamics to 
the analogues equations (36) and (37), respectively, used to derive "Hartree-Fock 
statics" . The former reduce to the latter if the dynamical aspect and the associated 
time evolution is neglected, that is in the limit that the nuclear and electronic 
momenta are absent or constant. Thus, the Car-Parrinello ansatz, namely Eq. (41) 
together with Eqs. (42)-(43), can also be viewed as a prescription to derive a new 
class of "dynamical ab initio methods" in very general terms. 



2-4-6 The Simulated Annealing and Optimization Viewpoints 

In the discussion given above, Car-Parrinello molecular dynamics was motivated 
by "combining" the positive features of both Ehrenfest and Born-Oppenheimer 
molecular dynamics as much as possible. Looked at from another side, the Car- 
Parrinello method can also be considered as an ingenious way to perform global 
optimizations (minimizations) of nonlinear functions, here (^ol^el^o), m a high 
dimensional parameter space including complicated constraints. The optimization 
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parameters are those used to represent the total wavefunction Hf in terms of simpler 
functions, for instance expansion coefficients of the orbitals in terms of Gaussians 
or plane waves, see e.g. Refs. 583 ' 375 ' 693 > 608 for applications of the same idea in 
other fields. 

Keeping the nuclei frozen for a moment, one could start this optimization pro- 
cedure from a "random wavefunction" which certainly does not minimize the elec- 
tronic energy. Thus, its fictitious kinetic energy is high, the electronic degrees of 
freedom are "hot". This energy, however, can be extracted from the system by 
systematically cooling it to lower and lower temperatures. This can be achieved 
in an elegant way by adding a non-conservative damping term to the electronic 
Car-Parrinello equation of motion Eq. (45) 

S 5 
tHi>i{t) = -T77 (*o|^o|*o) + Y77 {constraints} - j e fHipi , (61) 
dip* dip* 

where 7 > is a friction constant that governs the rate of energy dissipation 610 ; 
alternatively, dissipation can be enforced in a discrete fashion by reducing the ve- 
locities by multiplying them with a constant factor (1. Note that this deterministic 
and dynamical method is very similar in spirit to simulated annealing 332 invented 
in the framework of the stochastic Monte Carlo approach in the canonical ensemble. 
If the energy dissipation is done slowly, the wavefunction will find its way down to 
the minimum of the energy. At the end, an intricate global optimization has been 
performed! 

If the nuclei arc allowed to move according to Eq. (44) in the presence of an- 
other damping term a combined or simultaneous optimization of both electrons 
and nuclei can be achieved, which amounts to a "global geometry optimization" . 
This perspective is stressed in more detail in the review Ref. 223 and an imple- 
mentation of such ideas within the CADPAC quantum chemistry code is described in 
Ref. 692 . This operational mode of Car-Parrinello molecular dynamics is related to 
other optimization techniques where it is aimed to optimize simultaneously both the 
structure of the nuclear skeleton and the electronic structure. This is achieved by 
considering the nuclear coordinates and the expansion coefficients of the orbitals as 
variation parameters on the same footing 49 < 290 >608^ Car-Parrinello molecular 
dynamics is more than that because even if the nuclei continuously move according 
to Newtonian dynamics at finite temperature an initially optimized wavefunction 
will stay optimal along the nuclear trajectory. 

2-4-7 The Extended Lagrangian Viewpoint 

There is still another way to look at the Car-Parrinello method, namely in the 
light of so-called "extended Lagrangians" or "extended system dynamics" 14 , see 
e.g. Refs. 136 > 12 . 27 °! 585 > 217 fo r introductions. The basic idea is to couple additional 
degrees of freedom to the Lagrangian of interest, thereby "extending" it by increas- 
ing the dimensionality of phase space. These degrees of freedom are treated like 
classical particle coordinates, i.e. they arc in general characterized by "positions", 
"momenta", "masses", "interactions" and a "coupling term" to the particle's po- 
sitions and momenta. In order to distinguish them from the physical degrees of 
freedom, they are often called "fictitious degrees of freedom" . 
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The corresponding equations of motion follow from the Euler-Lagrange equa- 
tions and yield a microcanonical ensemble in the extended phase space where the 
Hamiltonian of the extended system is strictly conserved. In other words, the 
Hamiltonian of the physical (sub-) system is no more (strictly) conserved, and the 
produced ensemble is no more the microcanonical one. Any extended system dy- 
namics is constructed such that time-averages taken in that part of phase space that 
is associated to the physical degrees of freedom (obtained from a partial trace over 
the fictitious degrees of freedom) are physically meaningful. Of course, dynamics 
and thermodynamics of the system are affected by adding fictitious degrees of free- 
dom, the classic examples being temperature and pressure control by thermostats 
and barostats, see Sect. 4.2. 

In the case of Car-Parrincllo molecular dynamics, the basic Lagrangian for 
Newtonian dynamics of the nuclei is actually extended by classical fields {tpi(r)}, 
i.e. functions instead of coordinates, which represent the quantum wavefunction. 
Thus, vector products or absolute values have to be generalized to scalar products 
and norms of the fields. In addition, the "positions" of these fields {ipi} actually 
have a physical meaning, contrary to their momenta {ipi}- 

2.5 What about Hellmann-Feynman Forces ? 

An important ingredient in all dynamics methods is the efficient calculation of the 
forces acting on the nuclei, see Eqs. (30), (32), and (44). The straightforward 
numerical evaluation of the derivative 

Fj = -V/<*o|We|*o> (62) 

in terms of a finite-difference approximation of the total electronic energy is both 
too costly and too inaccurate for dynamical simulations. What happens if the gra- 
dients arc evaluated analytically? In addition to the derivative of the Hamiltonian 
itself 

V/ (* |We|*0> = (*o|Vj«e|*0> 

+ <Vj* |We|*0> + <*o|We|Vj* ) (63) 

there arc in general also contributions from variations of the wavefunction ~ V/ "Po- 
rn general means here that these contributions vanish exactly 

Ff FT = - (* |ViWe|* ) (64) 

if the wavefunction is an exact cigcnfunction (or stationary state wavefunction) of 
the particular Hamiltonian under consideration. This is the content of the often- 
cited Hcllmann-Feynman Theorem 295,186,368^ which is also valid for many varia- 
tional wavefunctions (e.g. the Hartree-Fock wavefunction) provided that complete 
basis sets are used. If this is not the case, which has to be assumed for numerical 
calculations, the additional terms have to be evaluated explicitly. 

In order to proceed a Slater determinant = det{V>i} of one-particle orbitals 
ipi, which themselves are expanded 

V^ = I>^(r;{R/}) (65) 
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in terms of a linear combination of basis functions {/„}, is used in conjunction with 
an effective one-particle Hamiltonian (such as e.g. in Hartree-Fock or Kohn-Sham 
theories). The basis functions might depend explicitly on the nuclear positions (in 
the case of basis functions with origin such as atom-centered orbitals), whereas the 
expansion coefficients always carry an implicit dependence. This means that from 
the outset two sorts of forces are expected 

V/V* = Yl ( V/C -) f^ r > { R '» + E c - ( v '/-( r ; { R '») ( 66 ) 

in addition to the Hcllmann- Fcynman force Eq. (64). 

Using such a linear expansion Eq. (65), the force contributions stemming from 
the nuclear gradients of the wavefunction in Eq. (63) can be disentangled into two 
terms. The first one is called "incomplete-basis-set correction" (IBS) in solid state 
theory 49 > 591 ! 180 and corresponds to the "wavefunction force" 494 or "Pulay force" in 
quantum chemistry 494 > 496 . It contains the nuclear gradients of the basis functions 

P? s = - ]T ((V,/, |K NSC - «| /„) + </„ K sc - c,| V/,)) (67) 

ivy 

and the (in practice non-self-consistent) effective one-particle Hamiltonian 49 » 591 . 
The second term leads to the so-called "non-self-consistency correction" (NSC) of 
the force 49 ' 591 

Ff c = -Jdr (V/n) (V SCF - V NSC ) (68) 

and is governed by the difference between the self-consistent ( "exact" ) potential or 
field V SCF and its non-self-consistent (or approximate) counterpart ]/ NSC associ- 
ated to Tt^ sc ; n(r) is the charge density. In summary, the total force needed in ab 
initio molecular dynamics simulations 

F, = Ff FT + Ff s + F? sc (69) 

comprises in general three qualitatively different terms; see the tutorial article 
Ref . 180 for a further discussion of core vs. valence states and the effect of pseudopo- 
tentials. Assuming that self-consistency is exactly satisfied (which is never going 
to be the case in numerical calculations), the force F^ sc vanishes and H^ CF has to 
be used to evaluate F j BS . The Pulay contribution vanishes in the limit of using a 
complete basis set (which is also not possible to achieve in actual calculations). 

The most obvious simplification arises if the wavefunction is expanded in terms 
of originless basis functions such as plane waves, see Eq. (100). In this case the Pu- 
lay force vanishes exactly, which applies of course to all ab initio molecular dynamics 
schemes (i.e. Ehrenfest, Born-Oppcnheimer, and Car-Parrincllo) using that par- 
ticular basis set. This statement is true for calculations where the number of plane 
waves is fixed. If the number of plane waves changes, such as in (constant pressure) 
calculations with varying cell volume / shape where the energy cutoff is strictly 
fixed instead, Pulay stress contributions crop up 219 > 245 > 660 > 211 > 202 ; se e Sect. 4.2. If 
basis sets with origin are used instead of plane waves Pulay forces arise always and 
have to be included explicitely in force calculations, see e.g. Refs. 75 > 370 > 371 for such 
methods. Another interesting simplification of the same origin is noted in passing: 
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there is no basis set superposition error (BSSE) 88 in plane wave-based electronic 
structure calculations. 

A non-obvious and more delicate term in the context of ab initio molecular 
dynamics is the one stemming from non-self-consistency Eq. (68) . This term van- 
ishes only if the wavefunction Wo is an eigenfunction of the Hamiltonian within the 
subspace spanned by the finite basis set used. This demands less than the Hcllmann- 
Feynman theorem where Wo has to be an exact eigenfunction of the Hamiltonian 
and a complete basis set has to be used in turn. In terms of electronic structure 
calculations complete self-consistency (within a given incomplete basis set) has to 
be reached in order that Fj sc vanishes. Thus, in numerical calculations the NSC 
term can be made arbitrarily small by optimizing the effective Hamiltonian and by 
determining its eigenfunctions to very high accuracy but it can never be suppressed 
completely. 

The crucial point is, however, that in Car-Parrinello as well as in Ehrcnfest 
molecular dynamics it is not the minimized expectation value of the electronic 
Hamiltonian, i.e. min* {(^ , o|W c |^'o)}, that yields the consistent forces. What is 
merely needed is to evaluate the expression (\f r o|W | , J r o) with the Hamiltonian and 
the associated wavefunction available at a certain time step, compare Eq. (32) to 
Eq. (44) or (30). In other words, it is not required (concerning the present discussion 
of the contributions to the force!) that the expectation value of the electronic 
Hamiltonian is actually completely minimized for the nuclear configuration at that 
time step. Whence, full self consistency is not required for this purpose in the case 
of Car-Parrinello (and Ehrenfcst) molecular dynamics. As a consequence, the non 
self-consistency correction to the force Fj sc Eq. (68) is irrelevant in Car-Parrinello 
(and Ehrenfest) simulations. 

In Born -Oppcnhcimcr molecular dynamics, on the other hand, the expectation 
value of the Hamiltonian has to be minimized for each nuclear configuration before 
taking the gradient to obtain the consistent force! In this scheme there is (inde- 
pendently from the issue of Pulay forces) always the non-vanishing contribution of 
the non-self-consistency force, which is unknown by its very definition (if it were 
know, the problem was solved, see Eq. (68)). It is noted in passing that there are 
estimation schemes available that correct approximately for this systematic error in 
Born-Oppenheimer dynamics and lead to significant time-savings, see e.g. Ref. 344 . 

Heuristically one could also argue that within Car-Parrinello dynamics the non- 
vanishing non-self-consistency force is kept under control or counterbalanced by 
the non-vanishing "mass times acceleration term" fiiipi(t) ss 0, which is small but 
not identical to zero and oscillatory. This is sufficient to keep the propagation sta- 
ble, whereas //ji/>i(f) = 0, i.e. an extremely tight minimization min* {(\I r o|?<o|\f r o)}, 
is required by its very definition in order to make the Born-Oppcnhcimcr approach 
stable, compare again Eq. (60) to Eq. (40). Thus, also from this perspective it 
becomes clear that the fictitious kinetic energy of the electrons and thus their ficti- 
tious temperature is a measure for the departure from the exact Born-Oppenheimer 
surface during Car-Parrinello dynamics. 

Finally, the present discussion shows that nowhere in these force derivations was 
made use of the Hellmann-Feynman theorem as is sometimes stated. Actually, it 
is known for a long time that this theorem is quite useless for numerical electronic 
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structure calculations, see e.g. Refs. 494 > 49 > 496 and references therein. Rather it 
turns out that in the case of Car-Parrinello calculations using a plane wave basis 
the resulting relation for the force, namely Eq. (64), looks like the one obtained by 
simply invoking the Hellmann-Feynman theorem at the outset. 

It is interesting to recall that the Hellmann-Feynman theorem as applied to a 
non-cigenfunction of a Hamiltonian yields only a first-order perturbative estimate 
of the exact force 295 > 368 . The same argument applies to ab initio molecular dy- 
namics calculations where possible force corrections according to Eqs. (67) and (68) 
are neglected without justification. Furthermore, such simulations can of course not 
strictly conserve the total Hamiltonian E cons Eq. (48). Finally, it should be stressed 
that possible contributions to the force in the nuclear equation of motion Eq. (44) 
due to position-dependent wavefunction constraints have to be evaluated following 
the same procedure. This leads to similar "correction terms" to the force, see e.g. 
Ref. 351 for such a case. 



2.6 Which Method to Choose ? 

Presumably the most important question for practical applications is which ab initio 
molecular dynamics method is the most efficient in terms of computer time given a 
specific problem. An a priori advantage of both the Ehrcnfest and Car-Parrincllo 
schemes over Born-Oppcnhcimcr molecular dynamics is that no diagonalization 
of the Hamiltonian (or the equivalent minimization of an energy functional) is 
necessary, except at the very first step in order to obtain the initial wavefunc- 
tion. The difference is, however, that the Ehrcnfest time evolution according to 
the time-dependent Schrodingcr equation Eq. (26) conforms to a unitary propaga- 
tion 3 4 ^366,342 

*(t + At) = exp[-iWe(t )At/fi]*(to) (70) 
y(t + mAt) = exp[-iH c (t + (m - I) At) At/h] 
x • • • 

x exp [-iH e (to + 2 At) At/h] 
x exp [-iH e {t + At) At/h] 
xexp[-iH c (t )At/h]y(t Q ) (71) 

t +t max 

dtH e (t) 



*(t + t max ) A *=°Texp 



to 



*(*o) (72) 



for infinitcsimally short times given by the time step At = i max /m; here T is the 
time -ordering operator and TL c {t) is the Hamiltonian (which is implicitly time 
dependent via the positions {Rj(t)}) evaluated at time t using e.g. split operator 
techniques 183 . Thus, the wavefunction W will conserve its norm and in particular 
orbitals used to expand it will stay orthonormal, see e.g. Ref. 617 . In Car-Parrinello 
molecular dynamics, on the contrary, the orthonormality has to be imposed brute 
force by Lagrange multipliers, which amounts to an additional orthogonalization 
at each molecular dynamics step. If this is not properly done, the orbitals will 
become non-orthogonal and the wavefunction unnormalizcd, see e.g. Sect. III.C.l 
in Ref. 472 . 
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But this theoretical disadvantage of Car-Parrinello vs. Ehrcnfcst dynamics is 
in reality more than compensated by the possibility to use a much larger time step 
in order to propagate the electronic (and thus nuclear) degrees of freedom in the 
former scheme. In both approaches, there is the time scale inherent to the nuclear 
motion T n and the one stemming from the electronic dynamics t c . The first one 
can be estimated by considering the highest phonon or vibrational frequency and 
amounts to the order of r n <~ 1CP 14 s (or 0.01 ps or 10 fs, assuming a maximum 
frequency of about 4000 cm -1 ). This time scale depends only on the physics of the 
problem under consideration and yields an upper limit for the timestep Ai max that 
can be used in order to integrate the equations of motion, e.g. At max w r n /10. 

The fasted electronic motion in Ehrenfest dynamics can be estimated within a 
plane wave expansion by ojf ~ E cu t, where i? cu t is the maximum kinetic energy 
included in the expansion. A realistic estimate for reasonable basis sets is ~ 
10~ 16 s, which leads to w r n /100. The analogues relation for Car-Parrinello 
dynamics reads however lu^ p <~ (E cnt / /i) 1 / 2 according to the analysis in Sect. 2.4, 
see Eq. (54). Thus, in addition to reducing by introducing a finite electron mass 
/i, the maximum electronic frequency increases much more slowly in Car-Parrinello 
than in Ehrenfest molecular dynamics with increasing basis set size. An estimate 
for the same basis set and a typical fictitious mass yields about t^ p ~ 10~ 15 s or 
T( p p w r n /10. According to this simple estimate, the time step can be about one 
order of magnitude larger if Car-Parrinello second-order fictitious-time electron 
dynamics is used instead of Ehrenfest first-order real-time electron dynamics. 

The time scale and thus time step problem inherent to Ehrenfest dynamics 
prompted some attempts to releave it. In Ref. 203 the equations of motion of 
electrons and nuclei were integrated using two different time steps, the one of the 
nuclei being 20-times as large as the electronic one. The powerful technology 
of multiple time step integration theory 636 > 639 could also be applied in order to 
ameliorate the time scale disparity 585 . A different approach borrowed from plasma 
simulations consists in decreasing the nuclear masses so that their time evolution 
is artificially speeded up 617 . As a result, the nuclear dynamics is fictitious (in the 
presence of real-time electron dynamics!) and has to be rescaled to the proper mass 
ratio after the simulation. 

In both Ehrenfest and Car-Parrinello schemes the explicitly treated electron 
dynamics limits the largest time step that can be used in order to integrate simul- 
taneously the coupled equations of motion for nuclei and electrons. This limitation 
does of course not exist in Born Oppenheimcr dynamics since there is no explicit 
electron dynamics so that the maximum time step is simply given by the one in- 
trinsic to nuclear motion, i.e. t®° « t„. This is formally an order of magnitude 
advantage with respect to Car-Parrinello dynamics. 

Do these back-of-the-envelope estimates have anything to do with reality? For- 
tunately, several state-of-the-art studies are reported in the literature for physi- 
cally similar systems where all three molecular dynamics schemes have been em- 
ployed. Ehrenfest simulations 553 > 203 of a dilute K^KCl)!-^ melt were performed 
using a time step of 0.012-0.024 fs. In comparison, a time step as large as 0.4 fs 
could be used to produce a stable Car-Parrinello simulation of electrons in liq- 
uid ammonia 155 ' 156 . Since the physics of these systems has a similar nature - 
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"unbound electrons" dissolved in liquid condensed matter (localizing as F'-centers, 
polarons, bipolarons, etc.) — the time step difference of about a factor of ten con- 
firms the crude estimate given above. In a Born-Oppcnhcimcr simulation 569 of 
again K^^KCl)!-^ but up to a higher concentration of unbound electrons the time 
step used was 0.5 fs. 

The time-scale advantage of Born-Oppcnhcimcr vs. Car-Parrincllo dynamics 
becomes more evident if the nuclear dynamics becomes fairly slow, such as in liquid 
sodium 343 or selenium 331 where a time step of 3 fs was used. This establishes 
the above-mentioned order of magnitude advantage of Born-Oppcnhcimcr vs. Car- 
Parrinello dynamics in advantageous cases. However, it has to be taken into account 
that in simulations 331 with such a large time step dynamical information is limited 
to about 10 THz, which corresponds to frequencies below roughly 500 cm -1 . In 
order to resolve vibrations in molecular systems with stiff covalent bonds the time 
step has to be decreased to less than a femtosecond (see the estimate given above) 
also in Born-Oppcnhcimcr dynamics. 

The comparison of the overall performance of Car-Parrincllo and Born- 
Oppcnhcimcr molecular dynamics in terms of computer time is a delicate issue. 
For instance it depends crucially on the choice made concerning the accuracy of 
the conservation of the energy i? con s as defined in Eq. (48). Thus, this issue is to 
some extend subject of "personal taste" as to what is considered to be a "suf- 
ficiently accurate" energy conservation. In addition, this comparison might to 
different conclusions as a function of system size. In order to nevertheless shed 
light on this point, microcanonical simulations of 8 silicon atoms were performed 
with various parameters using Car-Parrinello and Born-Oppenheimer molecular 
dynamics as implemented in the CPMD package 142 . This large-gap system was 
initially extremely well equilibrated and the runs were extended to 8 ps (and a 
few to 12 ps with no noticeable difference) at a temperature of about 360-370 K 
(with ±80 K root-mean-square fluctuations). The wavefunction was expanded up 
to E cut = 10 Ry at the T-point of a simple cubic supercell and LDA was used 
to describe the interactions. In both cases the velocity Verlet scheme was used to 
integrate the equations of motion, see Eqs. (231). It is noted in passing that also 
the velocity Verlet algorithm 638 allows for stable integration of the equations of 
motion contrary to the statements in Ref. 513 (see Sect. 3.4 and Figs. 4-5). 

In Car Parrincllo molecular dynamics two different time steps were used, 5 a.u. 
and 10 a.u. (corresponding to about 0.24 fs), in conjunction with a fictitious electron 
mass of fi — 400 a.u.; this mass parameter is certainly not optimized and thus 
the time step could be increased furthermore. Also the largest time step lead to 
perfect adiabaticity (similar to the one documented in Fig. 3), i.e. E p h ys Eq. (49) 
and T e Eq. (51) did not show a systematic drift relative to the energy scale set 
by the variations of V e Eq. (50). Within Born- Oppcnhcimcr molecular dynamics 
the minimization of the energy functional was done using the highly efficient DIIS 
(direct inversion in the iterative subspace) scheme using 10 "history vectors", see 
Sect. 3.6. In this case, the time step was either 10 a.u. or 100 a.u. and three 
convergence criteria were used; note that the large time step corresponding to 
2.4 fs is already at the limit to be used to investigate typical molecular systems 
(with frequencies up to 3-4000 cm -1 ). The convergence criterion is based on the 
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Figure 5. Conserved energy E conB defined in Eq. (48) from Car-Parrinello (CP) and Born- 
Oppcnhcimer (BO) molecular dynamics simulations of a model system for various time steps 
and convergence criteria using the CPMD package 142 ; see text for further details and Table 1 for 
the corresponding timings. Top: solid line: CP, 5 a.u.; open circles: CP, 10 a.u.; filled squares: 
BO, 10 a.u., 10~ 6 . Middle: open circles: CP, 10 a.u.; filled squares: BO, 10 a.u., 10~ 6 ; filled 
triangles: BO, 100 a.u., 10~ 6 ; open diamonds: BO, 100 a.u., 10~ 5 . Bottom: open circles: CP, 
10 a.u.; open diamonds: BO, 100 a.u., 10~ 5 ; dashed line: BO, 100 a.u., 10~ 4 . 



largest element of the wavefunction gradient which was required to be smaller than 
1CP 6 , 10~ 5 or 1CP 4 a.u.; note that the resulting energy convergence shows roughly 
a quadratic dependence on this criterion. 

The outcome of this comparison is shown in Fig. 5 in terms of the time evolution 
of the conserved energy E cons Eq. (48) on energy scales that cover more than three 
orders of magnitude in absolute accuracy. Within the present comparison ultimate 
energy stability was obtained using Car-Parrinello molecular dynamics with the 
shortest time step of 5 a.u., which conserves the energy of the total system to 
about 6xl0~ 8 a.u. per picosecond, see solid line in Fig. 5 (top). Increasing the 
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Table 1. Timings in CPU seconds and energy conservation in a.u. / ps for Car-Parrinello (CP) and 
Born-Oppcnhcimcr (BO) molecular dynamics simulations of a model system for 1 ps of trajectory 
on an IBM RS6000 / model 390 (Power2) workstation using the CPMD package 142 ; see Fig. 5 for 
corresponding energy plots. 



Method 


Time step (a.u.) 


Convergence (a.u.) 


Conservation (a.u./ps) 


Time (s) 


CP 


5 




6xl0" 8 


3230 


CP 


7 




lxl0~ 7 


2310 


CP 


10 




3xl0~ 7 


1610 


BO 


10 


10~ 6 


lxl0~ 6 


16590 


BO 


50 


10~ 6 


lxl0~ 6 


4130 


BO 


100 


10~ 6 


6xl0~ 6 


2250 


BO 


100 


10~ 5 


lxl0~ 5 


1660 


BO 


100 


io- 4 


lxlO" 3 


1060 



time step to 10 a.u. leads to an energy conservation of about 3xl0~ 7 a.u./ps and 
much larger energy fluctuations, see open circles in Fig. 5 (top). The computer time 
needed in order to generate one picosecond of Car-Parrinello trajectory increases - 
to a good approximation - linearly with the increasing time step, see Table 1 . The 
most stable Born-Oppenheimer run was performed with a time step of 10 a.u. and a 
convergence of 10~ 6 . This leads to an energy conservation of about 1 x 10~ 6 a.u./ps, 
see filled squares in Fig. 5(top). 

As the maximum time step in Born-Oppenheimer dynamics is only related 
to the time scale associated to nuclear motion it could be increased from 10 to 
100 a.u. while keeping the convergence at the same tight limit of 10~ 6 . This 
worsens the energy conservation slightly (to about 6xl0~ 6 a.u./ps), whereas the 
energy fluctuations increase dramatically, see filled triangles in Fig. 5(middle) and 
note the change of scale compared to Fig. 5(top). The overall gain is an acceleration 
of the Born-Oppenheimer simulation by a factor of about seven to eight, see Table 1. 
In the Born-Oppcnhcimcr scheme, the computer time needed for a fixed amount of 
simulated physical time decreases only sublincarly with increasing time step since 
the initial guess for the iterative minimization degrades in quality as the time step is 
made larger. Further savings of computer time can be easily achieved by decreasing 
the quality of the wavefunction convergence from 10~ 6 to 10~ 5 and finally to 10~ 4 , 
see Table 1. This is unfortunately tied to a significant decrease of the energy 
conservation from 6x 10~ 6 a.u./ps at 10~ 6 (filled triangles) to about 1 x 10 -3 a.u./ps 
at 10~ 4 (dashed line) using the same 100 a.u. time step, see Fig. 5(bottom) but 
note the change of scale compared to Fig. 5 (middle). 

In conclusion, Born-Oppcnhcimcr molecular dynamics can be made as fast 
as (or even faster than) Car-Parrinello molecular dynamics (as measured by the 
amount of CPU time spent per picosecond) at the expense of sacrificing accuracy 
in terms of energy conservation. In the "classical molecular dynamics community" 
there is a general consensus that this conservation law should be taken seriously 
being a measure of the numerical quality of the simulation. In the "quantum chem- 
istry and total energy communities" this issue is typically of less concern. There, it 
is rather the quality of the convergence of the wavefunction or energy (as achieved 
in every individual molecular dynamics step) that is believed to be crucial in order 
to gauge the quality of a particular simulation. 
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Finally, it is worth commenting in this particular section on a paper entitled 
"A comparison of Car-Parrinello and Born-Oppenhcimcr generalized valence bond 
molecular dynamics" 229 . In this paper one (computationally expensive) term in 
the nuclear equations of motion is neglected 648 > 405 . It is well known that using 
a basis set with origin, such as Gaussians /^(r;{R/}) centered at the nuclei, see 
Eq. (99), produces various Pulay forces, see Sect. 2.5. In particular a linear expan- 
sion Eq. (65) or (97) based on such orbitals introduces a position dependence into 
the orthogonality constraint 



that is hidden in the overlap matrix S'^ A1 ({R/}) which involves the basis functions. 
According to Eq. (44) this term produces a constraint force of the type 



in the correct Car-Parrinello equation of motion for the nuclei similar to the one 
contained in the electronic equation of motion Eq. (45). This term has to be 
included in order to yield exact Car-Parrinello trajectories and thus energy con- 
servation, see e.g. Eq. (37) in Ref. 351 for a similar situation. In the case of Born 
Oppcnhcimcr molecular dynamics, on the contrary, this term is always absent in the 
nuclear equation of motion, see Eq. (32). Thus, the particular implementation 229 
underlying the comparison between Car-Parrinello and Born Oppcnhcimcr molec- 
ular dynamics is an approximate one from the outset concerning the Car-Parrinello 
part; it can be argued that this was justified in the early papers 281 > 282 where the 
basic feasibility of both the Hartree Fock- and generalized valence bond-based Car- 
Parrinello molecular dynamics techniques was demonstrated 285 . Most importantly, 
this approximation implies that the energy E cons Eq. (48) cannot be rigorously con- 
served in this particular version of Car-Parrinello molecular dynamics. However, 
energy conservation of E cons was used in Ref. 229 to compare the efficiency and accu- 
racy of these two approaches to GVB ab initio molecular dynamics (using DIIS for 
the Born-Oppenheimer simulations as done in the above-given comparison). Thus, 
the final conclusion that for "... approaches that utilize non-space-fixed bases to 
describe the electronic wave function, Born-Oppenheimer AIMD is the method of 
choice, both in terms of accuracy and speed" cannot be drawn from this specific 
comparison for the reasons outlined above (independently of the particular basis 
set or electronic structure method used). 

The toy system investigated here (see Fig. 5 and Table 1), i.e. 8 silicon atoms in 
a periodic supcrccll, is for the purpose of comparing different approaches to ab initio 
molecular dynamics quite similar to the system used in Ref. 229 , i.e. clusters of 4 or 6 
sodium atoms (in addition, qualitatively identical results where reported in Sect. 4 
for silicon clusters). Thus, it is admissible to compare the energy conservations 
reported in Figs. 1 and 2 of Ref. 229 to the ones depicted here in Fig. 5 noting 
that the longest simulations reported in Ref. 229 reached only 1 ps. It should be 
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stressed that the energy conservation seen in Fig. 5 (top) is routinely achieved in 
Car-Parrincllo molecular dynamics simulations. 



2. 7 Electronic Structure Methods 

2. 7. 1 Introduction 

Up to this point, the electronic structure method to calculate the ab initio forces 
Vi^YHcY^) was not specified in detail. It is immediately clear that ab initio 
molecular dynamics is not tied to any particular approach, although very accu- 
rate techniques are of course prohibitively expensive. It is also evident that the 
strength or weakness of a particular ab initio molecular dynamics scheme is inti- 
mately connected to the strength or weakness of the chosen electronic structure 
method. Over the years a variety of different approaches such as density func- 
tional 108,679,35,472,343,36 ; H artree-Fock 365,254,i9i,379,28i,284,3i6,293 ; generalized va- 
lence bond (GVB) 282,283,228,229,230^ complete act i v e space SCF (CASSCF) 566 < 567 , 
full configuration interaction (FCI) 372 , scmicmpirical 669,671,91,190,114,666,280 or 
other approximate 473 < 454 > 551 > 4 55,i70,m,26 methods were combined with molecular 
dynamics, and this list is certainly incomplete. 

The focus of the present review clearly is Car-Parrinello molecular dynamics 
in conjunction with Hohenberg-Kohn-Sham density functional theory 301 > 338 . I n 
the following, only those parts of density functional theory are presented that im- 
pact directly on ab initio molecular dynamics. For a deeper presentation and in 
particular for a discussion of the assumptions and limitations of this approach 
(both conceptually and in practice) the reader is referred to the existing excellent 
literature 591 > 320 > 4 58,i68^ p or simplicity, the formulae are presented for the spin- 
unpolarized or restricted special case. 

Following the exposition of density functional theory, the fundamentals of 
Hartree-Fock theory, which is often considered to be the basis of quantum chem- 
istry, are introduced for the same special case. Finally, a glimpse is given at post 
Hartree-Fock methods. Again, an extensive text-book literature exists for these 
wavefunction-based approaches to electronic structure calculations 604 > 418 . The 
very useful connection between the density-based and wavefunction-based meth- 
ods goes back to Lowdin's work in the mid fifties and is e.g. worked out in Chapt. 2.5 
of Ref. 458 , where Hartree-Fock theory is formulated in density-matrix language. 

2.7.2 Density Functional Theory 

The total ground-state energy of the interacting system of electrons with classical 
nuclei fixed at positions {R/} can be obtained 

min{<* |«e|*o}} =minS KS [{0 4 }] 

*0 {ipi} 

as the minimum of the Kohn-Sham energy 301 > 338 

£ KS [{<M] = T„[{&}] + J dr V cxt (r) n(r) + \ J dr V H (r) n(r) + E KC [n] ,(75) 
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which is an explicit functional of the set of auxiliary functions {(pi(r)} that sat- 
isfy the orthonormality relation (<fii | (fij) = Sij. This is a dramatic simplification 
since the minimization with respect to all possible many-body wavefunctions {W} is 
replaced by a minimization with respect to a set of orthonormal one-particle func- 
tions, the Kohn Sham orbitals {4>i\. The associated electronic one-body density 
or charge density 

occ 

n(r) = ^/,|0 2 (r)| 2 (76) 

i 

is obtained from a single Slater determinant built from the occupied orbitals, where 
{fi} are integer occupation numbers. 

The first term in the Kohn-Sham functional Eq. (75) is the kinetic energy of a 
non-interacting reference system 



(77) 



consisting of the same number of electrons exposed to the same external potential 
as in the fully interacting system. The second term comes from the fixed external 
potential 



in which the electrons move, which comprises the Coulomb interactions between 
electrons and nuclei and in the definition used here also the internuclear Coulomb 
interactions; this term changes in the first place if core electrons are replaced by 
pseudopotentials, see Sect. 3.1.5 for further details. The third term is the Hartree 
energy, i.e. the classical electrostatic energy of two charge clouds which stem from 
the electronic density and is obtained from the Hartree potential 

V H (r) = J dr> , (79) 

which in turn is related to the density via 

V 2 V H (r) = -47rn(r) (80) 

Poisson's equation. The last contribution in the Kohn Sham functional, the 
exchange-correlation functional _E xc [n], is the most intricate contribution to the 
total electronic energy. The electronic exchange and correlation effects are lumped 
together and basically define this functional as the remainder between the exact 
energy and its Kohn-Sham decomposition in terms of the three previous contribu- 
tions. 

The minimum of the Kohn Sham functional is obtained by varying the energy 
functional Eq. (75) for a fixed number of electrons with respect to the density 
Eq. (76) or with respect to the orbitals subject to the orthonormality constraint, 
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see e.g. the discussion following Eq. (35) for a similar variational procedure. This 
leads to the Kohn-Sham equations 

\V 2 + Kxt(r) + V H (r) + Mr) = £ (81) 

j-iv 2 + T/ KS (r)| &(r) = ^A tJ ^(r) (82) 

ff^,(r)^A^(r), (83) 

which arc one-electron equations involving an effective one-particle Hamiltonian 
Hj? s with the local potential V KS . Note that H^ s nevertheless embodies the elec- 
tronic many-body effects by virtue of the exchange-correlation potential 

£>-ExcM 



Sn(r) 



= ^xc(r) • (84) 



A unitary transformation within the space of the occupied orbitals leads to the 
canonical form 

Hf s & = afr (85) 

of the Kohn-Sham equations, where {cj} are the eigenvalues. In conventional static 
density functional or "band structure" calculations this set of equations has to be 
solved self-consistently in order to yield the density, the orbitals and the Kohn 
Sham potential for the electronic ground state 487 . The corresponding total energy 
Eq. (75) can be written as 

£ KS = ]T - 1 / dr V u (r) n(r) + E xc [n] - J dr n(r ) , ( 86 ) 

where the sum over Kohn-Sham eigenvalues is the so-called "band-structure en- 
ergy" . 

Thus, Eqs. (81)-(83) together with Eqs. (39)-(40) define Born-Oppenhcimer 
molecular dynamics within Kohn Sham density functional theory, see e.g. 
Refs. 232,616,594,35,679,472,36,343,344 for guch implementations. The functional deriva- 
tive of the Kohn-Sham functional with respect to the orbitals, the Kohn-Sham 
force acting on the orbitals, can be expressed as 

^ = , (87) 

which makes clear the connection to Car-Parrincllo molecular dynamics, see 
Eq. (45). Thus, Eqs. (59)-(60) have to be solved with the effective one-particle 
Hamiltonian in the Kohn-Sham formulation Eqs. (81)-(83). In the case of Ehrcn- 
fest dynamics presented in Sect. 2.2, which will not be discussed in further detail 
at this stage, the Runge-Gross time-dependent generalization of density functional 
theory 258 has to be invoked instead, see e.g. Refs. 203 < 617 > 532 . 

Crucial to any application of density functional theory is the approximation of 
the unknown exchange and correlation functional. A discussion focussed on the 



35 



utilization of such functionals in the framework of ab initio molecular dynamics 
is for instance given in Ref. 588 . Those exchange-correlation functionals that will 
be considered in the implementation part, Sect. 3.3, belong to the class of the 
"Generalized Gradient Approximation" 

Eg GA [n] = Jdv n(r) sg GA (n(r); Vn(r)) , (88) 

where the unknown functional is approximated by an integral over a function that 
depends only on the density and its gradient at a given point in space, see Ref. 477 
and references therein. The combined exchange-correlation function is typically 
split up into two additive terms £ x and e c for exchange and correlation, respectively. 
In the simplest case it is the exchange and correlation energy density e^ A (n) of an 
interacting but homogeneous electron gas at the density given by the "local" density 
n(r) at space-point r in the inhomogencous system. This simple but astonishingly 
powerful approximation 320 is the famous local density approximation LDA 338 
(or local spin density LSD in the spin-polarized case ), and a host of different 
parameterizations exist in the literature 458 > 168 . The self-interaction correction 475 
SIC as applied to LDA was critically assessed for molecules in Ref. 240 with a 
disappointing outcome. 

A significant improvement of the accuracy was achieved by introducing the gra- 
dient of the density as indicated in Eq. (88) beyond the well-known straightforward 
gradient expansions. These so-called GGAs (also denoted as "gradient corrected" 
or "semilocal" functionals) extended the applicability of density functional calcula- 
tion to the realm of chemistry, see e.g. Refs. 476, 42, 362 ' 477 ' 478, 479 for a few "popular 
functionals" and Refs. 318 > 176 ,577,322 f or ex tensive tests on molecules, complexes, 
and solids, respectively. 

Another considerable advance was the successful introduction of "hybrid func- 
tionals" 43 ' 44 that include to some extent "exact exchange" 249 in addition to a 
standard GGA. Although such functionals can certainly be implemented within a 
plane wave approach 262 ' 128 ; they are prohibitively time-consuming as outlined at 
the end of Sect. 3.3. A more promising route in this respect are those function- 
als that include higher-order powers of the gradient (or the local kinetic energy 
density) in the sense of a generalized gradient expansion beyond the first term. 
Promising results could be achieved by including Laplacian or local kinetic energy 
terms 493 > 192 > 194 >662 > ^ ^ n j g s ^ a g e a sounc l judgment concerning their "prize / 
performance ratio" has to await further scrutinizing tests. The "optimized poten- 
tial method" (OPM) or "optimized effective potentials" (OEP) are another route 
to include "exact exchange" within density functional theory, see e.g. Sect. 13.6 
in Ref. 588 or Ref. 251 for overviews. Here, the exchange-correlation functional 
£° PM = i? xc [{0i}] depends on the individual orbitals instead of only on the den- 
sity or its derivatives. 

2.7.3 Hartree-Fock Theory 

Hartree-Fock theory is derived by invoking the variational principle in a restricted 
space of wavefunctions. The antisymmetric ground-state electronic wavefunction 
is approximated by a single Slater determinant = det{tpi} which is constructed 
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from a set of one-particle spin orbitals {ipi} required to be mutually orthonormal 
(ipi\ipj) — 5ij. The corresponding variational minimum of the total electronic 
energy H c denned in Eq. (2) 



~V 2 + V ext (r) 



1^2 J J drdr'^W^(r')|^^(r)^(r') 

\Y<j J drdr ' ^ (r) ^ (r ' } \^h\ ^ 



+ 



yields the lowest energy and the "best" wavefunction within a one-determinant 
ansatz; the external Coulomb potential V^xt was already defined in Eq. (78). Car- 
rying out the constraint minimization within this ansatz (see Eq. (36) in Sect. 2.3 
for a sketch) leads to 



[ 3 3) j 

-iv 2 + T/ HF (r)|^(r) = EA^W 



tf c HF ^(r) = E A ^( r ) 



(90) 

(91) 
(92) 



the Hartree-Fock integro-diffcrcntial equations. In analogy to the Kohn Sham 
equations Eqs. (81)-(83) these are effective one-particle equations that involve an 
effective one-particle Hamiltonian H^ F , the (Hartree-) Fock operator. The set of 
canonical orbitals 



is obtained similarly to Eq. (85). The Coulomb operator 



and the exchange operator 
/C, (r) Vi(r) 



(93) 



(94) 



(95) 



arc most easily defined via their action on a particular orbital Vi- It is found 
that upon acting on orbital ipi(r) the exchange operator for the j-th state "ex- 
changes" "0j( r ') ~~ * V'i( r ' / ) m the kernel as well as replaces ipi(r) — ► ?/>j(r) in its 
argument, compare to the Coulomb operator. Thus, /C is a non-local operator as 
its action on a function ipi at point r in space requires the evaluation and thus the 
knowledge of that function throughout all space by virtue of J dr' V>i( r ') • • • the 
required integration. In this sense the exchange operator does not possess a simple 
classical interpretation like the Coulomb operator C, which is the counterpart of 
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the Hartree potential Vh in Kohn-Sham theory. The exchange operator vanishes 
exactly if the antisymmetrization requirement of the wavefunction is relaxed, i.e. 
only the Coulomb contribution survives if a Hartree product is used to represent 
the wavefunction. 

The force acting on the orbitals is defined 

rpHF 

similarly to Eq. (87). At this stage, the various ab initio molecular dynamics 
schemes based on Hartrcc-Fock theory are defined, see Eqs. (39)-(40) for Born 
Oppcnhcimcr molecular dynamics and Eqs. (59)-(60) for Car-Parrinello molecu- 
lar dynamics. In the case of Ehrenfest molecular dynamics the time-dependent 
Hartree-Fock formalism 162 has to be invoked instead. 



2.7.4 Post Hartree-Fock Theories 

Although post Hartree-Fock methods have a very unfavorable scaling of the compu- 
tational cost as the number of electrons increases, a few case studies were performed 
with such correlated quantum chemistry techniques. For instance ab initio molec- 
ular dynamics was combined with GVB 282,283,228,229,230^ CASSCF 566,567^ as well 

as FCI 372 approaches, see also references therein. It is noted in passing that Car- 
Parrinello molecular dynamics can only be implemented straightforwardly if energy 
and wavefunction are "consistent". This is not the case in perturbation theories 
such as e.g. the widely used M0ller-Plesset approach 292 : within standard MP2 
the energy is correct to second order, whereas the wavefunction is the one given by 
the uncorrelated HF reference. As a result, the derivative of the MP 2 energy with 
respect to the wavefunction Eq. (96) does not yield the correct force on the HF 
wavefunction in the sense of fictitious dynamics. Such problems arc of course ab- 
sent from the Born-Oppcnhcimer approach to sample configuration space, see e.g. 
Ref. 328 > 317 > 33 for MP2, density functional, and multireference CI ab initio Monte 
Carlo schemes. 

It should be kept in mind that the rapidly growing workload of post HF calcu- 
lations, although extremely powerful in principle, limits the number of cxplicitcly 
treated electrons to only a few. The rapid development of correlated electronic 
structure methods that scale linearly with the number of electrons will certainly 
broaden the range of applicability of this class of techniques in the near future. 



2.8 Basis Sets 

2.8.1 Gaussians and Slater Functions 

Having selected a specific electronic structure method the next choice is related 
to which basis set to use in order to represent the orbitals tpi in terms of simple 
analytic functions f v with well-known properties. In general a linear combination 
of such basis functions 

^(r) = ^c w /„(r;{R / }) (97) 
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is used, which represents exactly any reasonable function in the limit of using a 
complete set of basis functions. In quantum chemistry, Slater-type basis functions 
(STOs) 

fl(r) = Nl r™»r™yr™* exp K m |r|] (98) 

with an exponentially decaying radial part and Gaussian-type basis functions 
(GTOs) 

/£(r) = Ng r™*r™*r?* exp [-a m r 2 ] (99) 

have received widespread use, see e.g. Rcf. 292 for a concise overview-type presen- 
tation. Here, N m , ( m and a m are constants that are typically kept fixed during 
a molecular electronic structure calculation so that only the orbital expansion co- 
efficients Ci v need to be optimized. In addition, fixed linear combinations of the 
above-given "primitive" basis functions can be used for a given angular momentum 
channel m, which defines the "contracted" basis sets. 

The Slater or Gaussian basis functions are in general centered at the positions of 
the nuclei, i.e. r — ► r R/ in Eq. (98)-(99), which leads to the linear combination 
of atomic orbitals (LCAO) ansatz to solve differential equations algebraically. Fur- 
thermore, their derivatives as well as the resulting matrix elements arc efficiently 
obtained by differentiation and integration in real-space. However, Pulay forces 
(see Sect. 2.5) will result for such basis functions that are fixed at atoms (or bonds) 
if the atoms are allowed to move, cither in geometry optimization or molecular 
dynamics schemes. This disadvantage can be circumvented by using freely floating 
Gaussians that are distributed in space 582 , which form an originlcss basis set since 
it is localized but not atom-fixed. 

2.8.2 Plane Waves 

A vastly different approach has its roots in solid-state theory. Here, the ubiquitous 
periodicity of the underlying lattice produces a periodic potential and thus imposes 
the same periodicity on the density (implying Bloch's Theorem, Born-von Karman 
periodic boundary conditions etc., see e.g. Chapt. 8 in Ref. 27 ). This heavily 
suggests to use plane waves as the generic basis set in order to expand the periodic 
part of the orbitals, see Sect. 3.1.2. Plane waves are defined as 

/P w (r)-7Vcxp[zGr] , (100) 

where the normalization is simply given by N = 1/ \/f2; fi is the volume of the 
periodic (super-) cell. Since plane waves form a complete and orthonormal set of 
functions they can be used to expand orbitals according to Eq. (97), where the 
labeling v is simply given by the vector G in reciprocal space / G-space (including 
only those G-vectors that satisfy the particular periodic boundary conditions). The 
total electronic energy is found to have a particularly simple form when expressed 
in plane waves 312 . 

It is important to observe that plane waves are originless functions, i.e. they 
do not depend on the positions of the nuclei {R/}- This implies that the Pulay 
forces Eq. (67) vanish exactly even within a finite basis (and using a fixed number 
of plane waves, see the discussion related to "Pulay stress" in Sect. 2.5), which 
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tremendously facilitates force calculations. This also implies that plane waves are 
a very unbiased basis set in that they are "delocalized" in space and do not "favor" 
certain atoms or regions over others, i.e. they can be considered as an ultimately 
"balanced basis set" in the language of quantum chemistry. Thus, the only way 
to improve the quality of the basis is to increase the "energy cutoff" E cut , i.e. to 
increase the largest |G|-vector that is included in the finite expansion Eq. (97). 
This blind approach is vastly different from the traditional procedures in quantum 
chemistry that are needed in order to produce reliable basis sets 292 . Another 
appealing feature is that derivatives in real-space are simply multiplications in G— 
space, and both spaces can be efficiently connected via Fast Fourier Transforms 
(FFTs). Thus, one can easily evaluate operators in that space in which they are 
diagonal, see for instance the flow charts in Fig. 6 or Fig. 7. 

According to the well-known "No Free Lunch Theorem" there cannot be only 
advantages connected to using plane waves. The first point is that the pseudopoten- 
tial approximation is intimately connected to using plane waves, why so? A plane 
wave basis is basically a lattice-symmetry-adapted three-dimensional Fourier de- 
composition of the orbitals. This means that increasingly large Fourier components 
are needed in order to resolve structures in real space on decreasingly small distance 
scales. But already orbitals of first row atoms feature quite strong and rapid oscilla- 
tions close to the nuclei due to the Pauli principle, which enforces a nodal structure 
onto the wavefunction by imposing orthogonality of the orbitals. However, most 
of chemistry is ruled by the valence electrons, whereas the core electrons are es- 
sentially inert. In practice, this means that the innermost electrons can be taken 
out of explicit calculations. Instead they are represented by a smooth and nodeless 
effective potential, the so-called pseudopotential 296,297,484,485,139^ gee £ or ms ^ ance 
Refs. 487,578,221 for revicws in tne contC xt of "solid state theory" and Refs. 145 ' 166 for 
pseudopotcntials as used in "quantum chemistry" . The resulting pseudo wavefunc- 
tion is made as smooth as possible close to the nuclear core region. This also means 
that properties that depend crucially on the wavefunction close to the core cannot 
be obtained straightforwardly from such calculations. In the field of plane wave 
calculations the introduction of "soft" norm-conserving ab initio pseudopotcntials 
was a breakthrough both conceptually 274 and in practice 28 . Another important 
contribution, especially for transition metals, was the introduction of the so-called 
ultrasoft pseudopotentials by Vanderbilt 661 . This approaches lead to the power- 
ful technique of plane wave-pseudopotential electronic structure calculations in the 
framework of density functional theory 312 > 487 . Within this particular framework 
the issue of pseudopotentials is elaborated in more detail in Sect. 3.1.5. 

Another severe shortcoming of plane waves is the backside of the medal of being 
an unbiased basis set: there is no way to shuffle more basis functions into regions in 
space where they are more needed than in other regions. This is particularly bad for 
systems with strong inhomogeneities. Such examples are all-electron calculations 
or the inclusion of semi-core states, a few heavy atoms in a sea of light atoms, and 
(semi-) finite systems such as surfaces or molecules with a large vacuum region in 
order to allow the long-range Coulomb interactions to decay. This is often referred 
to as the multiple length scale deficiency of plane wave calculations. 
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2.8.3 Generalized Plane Waves 



An extremely appealing and elegant generalization of the plane wave concept 263 > 264 
consists in denning them in curved £-space 

/g PW (0 = iVdct 1 / 2 Jcxp [iG r(0] (101) 
det J = 



where det J is the Jacobian of the transformation from Cartesian to curvilinear 
coordinates r — > £(r) with £ = (^ 1 ,^ 2 ,^ 3 ) and N = 1/y/U as for regular plane 
waves. These functions are orthonormal, form a complete basis set, can be used 
for k-point sampling after replacing G by G + k in Eq. (101), are originless (but 
nevertheless localized) so that Pulay forces are absent, can be manipulated via 
efficient FFT techniques, and reduce to standard plane waves in the special case of 
an Euclidean space £(r) = r. Thus, they can be used equally well like plane waves 
in linear expansions of the sort Eq. (65) underlying most of electronic structure 
calculations. The Jacobian of the transformation is related to the Riemannian 
metric tensor 

9iJ ~ 2^ t)r , dr j 
fe=l 

det J= det" 1/2 {#y} (102) 

which defines the metric of the £-space. The metric and thus the curvilinear co- 
ordinate system itself is considered as a variational parameter in the original fully 
adaptive-coordinate approach 263 . 264 ; see also Refs. 159 > 275 - 276 > 2 ?7,278^ xhus, a uni- 
form grid in curved Riemannian space is non-uniform or distorted when viewed in 
flat Euclidean space (where gij = Sij) such that the density of grid points (or the 
"local" cutoff energy of the expansion in terms of G-vectors) is highest in regions 
close to the nuclei and lowest in vacuum regions, sec Fig. 2 in Ref. 275 . 

Concerning actual calculations, this means that a lower number of generalized 
plane waves than standard plane waves arc needed in order to achieve a given ac- 
curacy 263 , see Fig. 1 in Ref. 275 . This allows even for all-electron approaches to 
electronic structure calculations where plane waves fail 431 > 497 . More recently, the 
distortion of the metric was frozen spherically around atoms by introducing defor- 
mation functions 265 ' 266 ; which leads to a concept closely connected to non uniform 
atom-centered meshes in real-space methods 431 , see below. In such non-fully- 
adaptive approaches using predefined coordinate transformations attention has to 
be given to Pulay force contributions which have to be evaluated explicitcly 265 > 431 . 



2.8.4 Wavelets 

Similar to using generalized plane waves is the idea to exploit the powerful 
multiscalc-properties of wavelets. Since this approach requires an extensive in- 
troductory discussion (see e.g. Ref. 242 for a gentle introduction) and since it seems 
still quite far from being used in large-scale electronic structure calculations the 
interested reader is referred to original papers 134 < 674 . 6 "<652,24i,25 anc j ^ ne g enera i 
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wavelet literature cited therein. Wavelet based methods allow intrinsically to ex- 
ploit multiple length scales without introducing Pulay forces and can be efficiently 
handled by fast wavelet transforms. In addition, they are also a powerful route to 
linear scaling or "order-TV" methods 453 ' 243 as first demonstrated in Ref. 241 with 
the calculation of the Hartree potential for an all electron uranium dimer. 



2.8.5 Mixed and Augmented Basis Sets 

Localized Gaussian basis functions on the one hand and plane waves on the other 
hand are certainly two extreme cases. There has been a tremendous effort to 
combine such localized and originless basis functions in order to exploit their mutual 
strengths. This resulted in a rich collection of mixed and augmented basis sets 
with very specific implementation requirements. This topic will not be covered 
here and the interested reader is referred to Refs. 75 > 654 . 498 > 370 > 371 anc j references 
given therein for some recent implementations used in conjunction with ah initio 
molecular dynamics. 



2.8.6 Wannier Functions 

An alternative to the plane wave basis set in the framework of periodic calculations 
in solid-state theory are Wannier functions, see for instance Sect. 10 in Ref. 27 . 
These functions are formally obtained from a unitary transformation of the Bloch 
orbitals Eq. (114) and have the advantage that they can be exponentially localized 
under certain circumstances. The so-called maximally localized generalized Wan- 
nier functions 413 are the periodic analogues of Boys' localized orbitals defined for 
isolated systems. Recently the usefulness of Wannier functions for numerical pur- 
poses was advocated by several groups, see Refs. 339 > 184 < 413 > 10 anc i references given 
therein. 



2.8.7 Real Space Grids 



A quite different approach is to leave conventional basis set approaches altogether 
and to resort to real-space methods where continuous space is replaced by a discrete 
space r — > r p . This entails that the derivative operator or the entire energy ex- 
pression has to be discretized in some way. The high-order central finite difference 
approach leads to the expression 



-jV Mr) = -J 



EN 



-AT ^n, *Pi( r Px i r p v + n y h, r p ^ ) 



+ T,Z=-N C n* i>i{r Px ,r Py ,r Pz +n z h) 



+ 0(h 2N+2 ) (103) 



for the Laplacian which is correct up to the order h 2N+2 . Here, h is the uniform 
grid spacing and {C n } are known expansion coefficients that depend on the selected 
order 130 . Within this scheme, not only the grid spacing h but also the order are 
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disposable parameters that can be optimized for a particular calculation. Note that 
the discretization points in continuous space can also be considered to constitute a 
sort of "finite basis set" - despite different statements in the literature - and that 
the "infinite basis set limit" is reached as h — > for N fixed. A variation on the 
theme are Mehrstellen schemes where the discretization of the entire differential 
equation and not only of the derivative operator is optimized 89 . 

The first real-space approach devised for ah initio molecular dynamics was based 
on the lowest-order finite-difference approximation in conjunction with a equally- 
spaced cubic mesh in real space 109 . A variety of other implementations of more 
sophisticated real-space methods followed and include e.g. non-uniform meshes, 
multigrid acceleration, different discretization techniques, and finite-element meth- 
ods 686,61,39,130,131,632,633,431,634 Among the chief advantages of the real-space 
methods is that linear scaling approaches 453 > 243 can be implemented in a natural 
way and that the multiple-length scale problem can be coped with by adapting the 
grid. However, the extension to such non uniform meshes induces the (in)famous 
Pulay forces (see Sect. 2.5) if the mesh moves as the nuclei move. 

3 Basic Techniques: Implementation within the CPMD Code 

3.1 Introduction and Basic Definitions 

This section discusses the implementation of the plane wave-pseudopotential molec- 
ular dynamics method within the CPMD computer code 142 . It concentrates on the 
basics leaving advanced methods to later chapters. In addition all formulas are for 
the non-spin polarized case. This allows to show the essential features of a plane 
wave code as well as the reasons for its high performance in detail. The imple- 
mentation of other versions of the presented algorithms and of the more advanced 
techniques in Sect. 4 is in most cases very similar. 

There are many reviews on the pseudopotential plane wave method alone or in 
connection with the Car-Parrincllo algorithm. Older articles 312 ' 157 > 487 > 591 as well 
as the book by Singh 578 concentrate on the electronic structure part. Other re- 
views 513 < 472 >223,224 p resen t the plane wave method in connection with the molecular 
dynamics technique. 

3.1.1 Unit Cell and Plane Wave Basis 

The unit cell of a periodically repeated system is defined by the Bravais lattice 
vectors ai, a 2 , and a 3 . The Bravais vectors can be combined into a three by three 
matrix h = [ai , a 2 , a 3 ] 459 . The volume f2 of the cell is calculated as the determinant 
of h 

S! = dcth . (104) 
Further, scaled coordinates s are introduced that are related to r via h 

r = hs . (105) 
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Distances in scaled coordinates are related to distances in real coordinates by the 
metric tensor Q = h*h 

( r< - rjf = ( 8j - stfGisi - Bj ) . (106) 

Periodic boundary conditions can be enforced by using 

r pbc = r-h[h _1 r] NINT , (107) 

where [•••]nint denotes the nearest integer value. The coordinates r p b c will be 
always within the box centered around the origin of the coordinate system. Recip- 
rocal lattice vectors bj are defined as 

bj • a,- = 27r<5j,- (108) 

and can also be arranged to a three by three matrix 

[bi,b 2 ,b 3 ] -2 7 r(h t )- 1 . (109) 

Plane waves build a complete and orthonormal basis with the above periodicity 
(see also the section on plane waves in Sect. 2.8) 

/P w (r) = -^=cxp[»G.r] = -i=exp[2^g.s] , (110) 

with the reciprocal space vectors 

G = 2 7 r(h t )- 1 g , (HI) 

where g = k] is a triple of integer values. A periodic function can be expanded 
in this basis 

^(r)=^(r + L) = -^£>(G)exp[iG-r] , (112) 

where tp(r) and -0(G) are related by a three-dimensional Fourier transform. The 
direct lattice vectors L connect equivalent points in different cells. 

3.1.2 Plane Wave Expansions 

The Kohn Sham potential (see Eq. (82)) of a periodic system exhibits the same 
periodicity as the direct lattice 

V KS (r) = V KS (r + L) , (113) 

and the Kohn-Sham orbitals can be written in Bloch form (see e.g. Ref. 27 ) 

*(r) = *j(r,k) =exp[ik-r] Uj(r,k) , (114) 

where k is a vector in the first Brillouin zone. The functions Wj(r,k) have the 
periodicity of the direct lattice 

Ui(r,k) = u,(r + L,k) . (115) 
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The index i runs over all states and the states have an occupation /j(k) associated 
with them. The periodic functions u,(r, k) are now expanded in the plane wave 
basis 

Ul (r,k) = -^^ Cj (G,k)exp[jG-r] , (116) 

G 

and the Kohn-Sham orbitals are 

<^(r,k) = -i=^ Ci (G,k)exp[i(G + k).r] , (117) 

where Cj(G,k) are complex numbers. With this expansion the density can also be 
expanded into a plane wave basis 

= f^fiM E c*(G',k) Cl (G,k)cxp[z(G + k)-r] (118) 

i J G,G' 

= ^n(G)cxp[iG t] , (119) 



where the sum over G vectors in Eq. (119) expands over double the range given 
by the wavefunction expansion. This is one of the main advantages of the plane 
wave basis. Whereas for atomic orbital basis sets the number of functions needed 
to describe the density grows quadratically with the size of the system, there is 
only a linear dependence for plane waves. 

3.1.3 K-Points and Cutoffs 

In actual calculations the infinite sums over G vectors and cells has to be truncated. 
Furthermore, we have to approximate the integral over the Brillouin zone by a finite 
sum over special k-points 

fdk -^J2w k , (120) 

•* k 

where Wk are the weights of the integration points. Schemes on how to choose the 
integration points efficiently are available in the literature 30 > 123 > 435 where also an 
overview 179 on the use of k-points in the calculation of the electronic structure of 
solids can be found. 

The truncation of the plane wave basis rests on the fact that the Kohn-Sham 
potential V KS (G) converges rapidly with increasing modulus of G. For this reason, 
at each k-point, only G vectors with a kinetic energy lower than a given maximum 
cutoff 

X - |k + G| 2 <E cut (121) 

are included in the basis. With this choice of the basis the precision of the calcu- 
lation within the approximations of density functional theory is controlled by one 
parameter E cut only. 
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The number of plane waves for a given cutoff depends on the unit cell and the 
k-point. An estimate for the size of the basis at the center of the Brillouin zone is 

AW = , (122) 

where E cut is in Hartree units. The basis set needed to describe the density cal- 
culated from the Kohn-Sham orbitals has a corresponding cutoff that is four times 
the cutoff of the orbitals. The number of plane waves needed at a given density 
cutoff is therefore eight times the number of plane waves needed for the orbitals. 

It is a common approximation in density functional theory calculations 536 > 169 
to use approximate electronic densities. Instead of using the full description, the 
density is expanded in an auxiliary basis. An incomplete plane wave basis can be 
considered as an auxiliary basis with special properties 371 . Because of the filter 
property of plane waves the new density is an optimal approximation to the true 
density. No additional difficulties in calculations of the energy or forces appear. 
The only point to control is, if the accuracy of the calculation is still sufficient. 

Finally, sums over all unit cells in real space have to be truncated. The only 
term in the final energy expression with such a sum is the real space part of the 
Ewald sum (see Sect. 3.2). This term is not a major contribution to the workload 
in a density functional calculation, that is the cutoff can be set rather generously. 



3.1.4 Real Space Grid 

A function given as a finite linear combination of plane waves can also be defined 
as a set of functional values on a equally spaced grid in real space. The sampling 
theorem (see e.g. Ref. 492 ) gives the maximal grid spacing that still allows to hold 
the same information as the expansion coefficients of the plane waves. The real 
space sampling points R are defined 

R= hNq , (123) 

where N is a diagonal matrix with the entries 1/N S and q is a vector of integers 
ranging from to N s — 1 (s = x, y, z). To fulfill the sampling theorem N s has to 
be bigger than 2max(g s ) + 1. To be able to use fast Fourier techniques, N s must 
be decomposable into small prime numbers (typically 2, 3, and 5). In applications 
the smallest number N s that fulfills the above requirements is chosen. 
A periodic function can be calculated at the real space grid points 

/(R) = £/(G)exp[iG-R] (124) 

G 

= ^/(G)exp[2 7 r»((h t )- 1 g) • (hNq)] (125) 
g 

= £ /(G) exp 
g 

The function /(G) is zero outside the cutoff region and the sum over g can be 
extended over all indices in the cube — g™ ax . . . g™ ax . The functions /(R) and 
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/(G) arc related by three-dimensional Fourier transforms 



/(R) = inv_FT [/(G)] 
/(G) = fw_FT [/(R)] . 

The Fourier transforms are defined by 

Wx-l N y -l at z _i 



(127) 
(128) 
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where the appropriate mappings of q and g to the indices 
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(130) 



(131) 
(132) 
(133) 



have to be used. From Eqs. (129) and (130) it can be seen, that the calculation 
of the three-dimensional Fourier transforms can be performed by a series of one 
dimensional Fourier transforms. The number of transforms in each direction is 
N x N y , N X N Z , and N y N z respectively. Assuming that the one-dimensional trans- 
forms are performed within the fast Fourier transform framework, the number of 
operations per transform of length n is approximately 5nlogn. This leads to an 
estimate for the number of operations for the full three-dimensional transform of 
BNlogN, where N = N x N y N z . 



3.1.5 Pseudopotentials 

In order to minimize the size of the plane wave basis necessary for the calculation, 
core electrons are replaced by pseudopotentials. The pscudopotcntial approxima- 
tion in the realm of solid-state theory goes back to the work on orthogonalizcd 
plane waves 298 and core state projector methods 485 . Empirical pseudopotentials 
were used in plane wave calculations 294 ' 703 but new developments have consid- 
erably increased efficiency and reliability of the method. Pseudopotential are re- 
quired to correctly represent the long range interactions of the core and to produce 
pseudo-wavefunction solutions that approach the full wavefunction outside a core 
radius r c . Inside this radius the pseudopotential and the wavefunction should be as 
smooth as possible, in order to allow for a small plane wave cutoff. For the pseudo- 
wavefunction this requires that the nodal structure of the valence wavefunctions 
is replaced by a smooth function. In addition it is desired that a pseudopotential 
is transferable 238 > 197 ; this means that one and the same pseudopotential can be 
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used in calculations of different chemical environment resulting in calculations with 
comparable accuracy. 

A first major step to achieve all this conflicting goals was the introduction of 
"norm-conservation" 622 ' 593 . Norm-conserving pseudopotentials have to be angular 
momentum dependent. In their most general form they are semi-local 



V PF (r,r') = J2Ylm(r)V l (r)S T yY lm (r') , (134) 

lm 



where Y\_ m are spherical harmonics. A minimal set of requirements and a construc- 
tion scheme for soft, semi-local pseudopotentials were developed 274 > 28 . Since then 
many variations of the original method have been proposed, concentrating cither 
on an improvement in softness or in transferability. Analytic representations of the 
core part of the potential 326 . 626, 627,509 were used. Extended norm-conservation 564 
was introduced to enhance transferability and new concepts to increase the soft- 
ness were presented 659 > 509 > 369 . More information on pseudopotentials and their 
construction can be found in recent review articles 487 < 578 > 221 . 

Originally generated in a semi-local form, most applications use the fully separa- 
ble form. Pseudopotentials can be transformed to the separable form using atomic 
wavefunctions 335 > 73 > 659 . Recently 239 > 288 a new type of separable, norm-conserving 
pseudopotentials was introduced. Local and non-local parts of these pseudopoten- 
tials have a simple analytic form and only a few parameters are needed to charac- 
terize the potential. These parameters are globally optimized in order to reproduce 
many properties of atoms and ensure a good transferability. 

A separable non-local pseudopotential can be put into general form (this in- 
cludes all the above mentioned types) 

t/ pp (r,r') = (V cole (r) + AV local (r)) S T , r , + Y,Pk(r)hklPl(r') • (135) 

k,l 

The local part has been split into a core ( ~ 1 jr for r — > oo ) and a short-ranged 
local part in order to facilitate the derivation of the final energy formula. The 
actual form of the core potential will be defined later. The local potential AV[ oca i 
and the projectors Pk are atom-centered functions of the form 

< P (t) = <p(\t-B. i \) Y lm (e,<j>) , (136) 

that can be expanded in plane waves 

p(r) = ]T <p{G) exp[iG • r] S/(G) Y lm (6, $) , (137) 

G 

R/ denote atomic positions and the so-called structure factors Si are defined as 

S/(G) = expHG • R 7 ] . (138) 
The functions <p(G) are calculated from ip(r) by a Bessel transform 

/■OO 

<p(G) = 4n (-i) 1 drr 2 p(r)ji(Gr) , (139) 
Jo 
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where j; are spherical Bessel functions of the first kind. The local pseudopotcntial 
and the projectors of the nonlocal part in Fourier space are given by 



AMocal(G) 
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drr 2 A Viocai (r)jo (Gr) 



(140) 



P k (G) = -j= (-i) 1 j drr'P^MGr) Y lm (6,4>) , (141) 
where Im are angular momentum quantum numbers associated with projector a. 

3.2 Electrostatic Energy 
3.2.1 General Concepts 

The electrostatic energy of a system of nuclear charges Zj at positions R/ and 
an electronic charge distribution n(r) consists of three parts: the Hartree energy 
of the electrons, the interaction energy of the electrons with the nuclei and the 
intcrnuclcar interactions 



drdv 1 



./ n(r)n(r') 
Ir-r'l 



i J i^j 1 J 1 



The Ewald method (see e.g. Ref. 12 ) can be used to avoid singularities in the 
individual terms when the system size is infinite. In order to achieve this a Gaussian 
core charge distribution associated with each nuclei is defined 



n,(r) = o 

(Rj) 3 



7r 3/2 exp 



R; 



(143) 



It is convenient at this point to make use of the arbitrariness in the definition of the 
core potential and define it to be the potential of the Gaussian charge distribution 
of Eq. (143) 



V c { vc (r) = J dr' 



Zi 



Ir-Rj 



-erf 



R/ 



(144) 



where erf is the error function. The interaction energy of this Gaussian charge 
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distributions is now added and subtracted from the total electrostatic energy 

, n(r)n(r') 



Ees 



dvdv' 



drdr' 
dvdv 



r - r'| 
, n c (r)n c (r') 
r-r'| 
, rc c (r)n(r') 



1 \ ~> ZjZj 
2£-, |Rj-Rj| 

, ra c (r)n c (r') 



drdr' 



(145) 



where n c (r) — J2i n l( r )- The first three terms can be combined to the electrostatic 
energy of a total charge distribution n to t(r) = n(r) + n c (r). The remaining terms 
are rewritten as a double sum over nuclei and a sum over self-energy terms of the 
Gaussian charge distributions 



Ees= 2 



drdr 



, n t ot(r)ntot(r') 



1 ZjZj 

x > tt; — — rCrfc 



IRj -R./| 



E 



2ir Rj 



(146) 



where erfc denotes the complementary error function. 



3.2.2 Periodic Systems 

For a periodically repeated system the total energy per unit cell is derived from 
the above expression by using the solution to Poisson's equation in Fourier space 
for the first term and make use of the quick convergence of the second term in real 
space. The total charge is expanded in plane waves with expansion coefficients 



n tot (G) = n(G) + ]T ^(G)Sj(G) 
i 



n(G) 



exp 



\g 2 r? 



Si(G) 



This leads to the electrostatic energy for a periodic system 

Ees — 2?r ^ ! to jj; 2 - + £ ovrl - £ se if 
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and 

Here, the sums expand over all atoms in the simulation cell, all direct lattice vectors 
L, and the prime in the first sum indicates that I(J is imposed for L = 0. 

3.2.3 Cluster Boundary Conditions 

The possibility to use fast Fourier transforms to calculate the electrostatic energy 
is one of the reasons for the high performance of plane wave calculations. How- 
ever, plane wave based calculations imply periodic boundary conditions. This is 
appropriate for crystal calculations but very unnatural for molecule or slab calcu- 
lations. For neutral systems this problem is circumvented by use of the supercell 
method. Namely, the molecule is periodically repeated but the distance between 
each molecule and its periodic images is so large that their interaction is negligible. 
This procedure is somewhat wasteful but can lead to satisfactory results. 

Handling charged molecular systems is, however, considerably more difficult, 
due to the long range Coulomb forces. A charged periodic system has infinite 
energy and the interaction between images cannot really be completely eliminated. 
In order to circumvent this problem several solutions have been proposed. The 
simplest fix-up is to add to the system a neutralizing background charge. This 
is achieved trivially as the G — term in Eq. (149) is already eliminated. This 
leads to finite energies but does not eliminate the interaction between the images 
and makes the calculation of absolute energies difficult. Other solutions involve 
performing a set of different calculations on the system such that extrapolation to 
the limit of infinitely separated images is possible. This procedure is lengthy and 
one cannot use it easily in molecular dynamics applications. It has been shown, 
that it is possible to estimate the correction to the total energy for the removal 
of the image charges 378 . Still it seems not easy to incorporate this scheme into 
the frameworks of molecular dynamics. Another method 60 > 702 > 361 works with the 
separation of the long and short range parts of the Coulomb forces. In this method 
the low-order multipole moments of the charge distribution are separated out and 
handled analytically. This method was used in the context of coupling ab initio 
and classical molecular dynamics 76 . 

The long-range forces in Eq. (146) are contained in the first term. This term 
can be written 

\ jjdrdr' > = \ J d rV H (r)n tot (r) , (152) 

where the electrostatic potential Vh(i") is the solution of Poisson's equation (see 
Eq. (80)). There are two approaches to solve Poisson's equation subject to the 
boundary conditions Va(r) — > for r — > oo implemented in CPMD. Both of them 
rely on fast Fourier transforms, thus keeping the same framework as for the periodic 
case. 

The first method is due to Hockney 300 and was first applied to density functional 
plane wave calculations in Ref. 36 . In the following outline, for the sake of simplicity, 
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a one-dimensional case is presented. The charge density is assumed to be non-zero 
only within an interval L and sampled on N equidistant points. These points are 
denoted by x p . The potential can then be written 

L °° 

Vk(ttp) = — ^2 G ( x p - v) n 0v) (153) 
p'—— 00 

L N 

= XI g ( x p ~ vMv) ( 154 ) 

p'=0 

for p = 0, 1, 2, . . . N, where G(x p — x p >) is the corresponding Green's function. In 
Hockney's algorithm this equation is replaced by the cyclic convolution 

L 2N+1 

Vh(x p ) = — ^2 G(x p - x p >)n(x p >) (155) 

p'=0 

where p = 0,1,2, ...2N+1, and 

nix ) - <p<iV 

n ^P>-\0 N<p<2N+l (i5bj 

G(x p ) = G(x p ) - (N + 1) < p < N (157) 
h(xp) = n(x p + L) (158) 
G(xp) = G(x p + L) (159) 

The solution Vh(x p ) can be obtained by a series of fast Fourier transforms and has 
the desired property 

Vk(x p ) = y H (x p ) for < p < N . (160) 

To remove the singularity of the Green's function at x = 0, G(a;) is modified for 
small x and the error is corrected by using the identity 



G(x) = -erf 
x 



X 




X 




+ — crfc 






X 





(161) 



where r c is chosen such, that the short-ranged part can be accurately described by 
a plane wave expansion with the density cutoff. In an optimized implementation 
Hockney's method requires the double amount of memory and two additional fast 
Fourier transforms on the box of double size (see Fig. 6 for a flow chart). Hockney's 
method can be generalized to systems with periodicity in one (wires) and two (slabs) 
dimensions. It was pointed out 173 that Hockney's method gives the exact solution 
to Poisson's equation for isolated systems if the boundary condition (zero density 
at the edges of the box) are fulfilled. 

A different, fully reciprocal space based method, that can be seen as an approx- 
imation to Hockney's method, was recently proposed 393 . The final expression for 
the Hartree energy is also based on the splitting of the Green's function in Eq. (161) 

£ ES = 2^^ V^ T (G)nt ot (G) + E ovll - E scli . (162) 
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Figure 6. Flow chart for the calculation of long-ranged part of the electrostatic energy using the 
method by Hockney 300 . The part inside the dashed box is calculated most efficiently with the 
procedure outlined by Eastwood and Brownrigg 173 . 



The potential function is calculated from two parts, 

< T (G) = Vk(G) + Vk(G) , (163) 
where Vh(G) is the analytic part, calculated from a Fourier transform of erfc 
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Vk(G) = — ( 1 - exp 



G 2 rl 



n(G) 



(164) 



and Vh(G) is calculated from a discrete Fourier transform of the Green's function 
on an appropriate grid. The calculation of the Green's function can be done at the 
beginning of the calculation and has not to be repeated again. It is reported 393 
that a cutoff of ten to twenty percent higher than the one employed for the charge 
density gives converged results. The same technique can also be applied for systems 
that are periodic in one and two dimensions 394 . 

If the boundary conditions are appropriately chosen, the discrete Fourier trans- 
forms for the calculation of Vh(G) can be performed analytically 437 . This is 
possible for the limiting case where r c = and the boundary conditions arc on a 
sphere of radius R for the cluster. For a one-dimensional system we choose a torus 
of radius R and for the two-dimensional system a slab of thickness Z . The clcctro- 

1/2 

static potential for these systems are listed in Table 2, where G xy = [g* + <?y] 
and J n and K „ are the Bessel functions of the first and second kind of integer order 
n. 

Hockney's method requires a computational box such that the charge density is 
negligible at the edges. This is equivalent to the supercell approach 510 . Practical 
experience tells that a minimum distance of about 3 A of all atoms to the edges of 
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Table 2. Fourier space formulas for the Hartree energy, see text for definitions. 
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the box is sufficient for most systems. The Green's function is then applied to the 
charge density in a box double this size. The Green's function has to be calculated 
only once at the beginning of the calculation. The other methods presented in this 
chapter require a computational box of double the size of the Hockncy method as 
they are applying the artificially periodic Green's function within the computational 
box. This can only be equivalent to the exact Hockney method if the box is enlarged 
to double the size. In plane wave calculations computational costs grow linearly 
with the volume of the box. Therefore Hockney's method will prevail over the others 
in accuracy, speed, and memory requirements in the limit of large systems. The 
direct Fourier space methods have advantages through their easy implementation 
and for small systems, if not full accuracy is required, i.e. if they are used with 
smaller computational boxes. In addition, they can be of great use in calculations 
with classical potentials. 



3.3 Exchange and Correlation Energy 

Exchange and correlation functionals implemented in the CPMD code are of the local 
type with gradient corrections. These type of functionals can be written as (see 
also Eqs. (88) and (84)) 



£ xc = / dre xc (n,Vn)n(r) = VL ^ e xc (G)n*(G) 



with the corresponding potential 

OF- 

On ^ Or 



T/ , s dF xc ^ 
^xc(r) = ^--^ 



0{0 s n) 



(165) 



(166) 



where _F XC = e xc (n, Vn)n and d s n is the s-componcnt of the density gradient. 

Exchange and correlation functionals have complicated analytical forms that 
give rise to high frequency components in e xc (G). Although these high frequency 
components do not enter the sum in Eq. (165) due to the filter effect of the density, 
they affect the calculation of e xc . As the functionals are only local in real space, not 
in Fourier space, they have to be evaluated on a real space grid. The function e xc (G) 
can then be calculated by a Fourier transform. Therefore the exact calculation of 
E xc would require a grid with infinite resolution. However, the high frequency 
components are usually very small and even a moderate grid gives accurate results. 
The use of a finite grid results in an effective redefinition of the exchange and 
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correlation energy 

Exc = ~NWW £ £xc(n ' V ")( R )«( R ) = n E MG)n(G) , (167) 

x y z R G 

where e xc (G) is the finite Fourier transform of e xc (R). This definition of £ xc 
allows the calculation of all gradients analytically. In most applications the real 
space grid used in the calculation of the density and the potentials is also used 
for the exchange and correlation energy. Grids with higher resolution can be used 
easily. The density is calculated on the new grid by use of Fourier transforms and 
the resulting potential is transfered back to the original grid. With this procedure 
the different grids do not have to be commensurate. 

The above redefinition has an undesired side effect. The new exchange and 
correlation energy is no longer translationally invariant. Only translations by a 
multiple of the grid spacing do not change the total energy. This introduces a 
small modulation of the energy hyper surface 685 , known as " ripples" . Highly 
accurate optimizations of structures and the calculation of harmonic frequencies 
can be affected by the ripples. Using a denser grid for the calculation of E xc is the 
only solution to avoid these problems. 

The calculation of a gradient corrected functional within the plane wave frame- 
work can be conducted using Fourier transforms 685 . The flowchart of the calcula- 
tion is presented in Fig. 7. With the use of Fourier transforms the calculation of 
second derivatives of the charge density is avoided, leading to a numerically stable 
algorithm. To this end the identity 

dF xc dF xc d s n 

d(d s n) ~ d\Wn\ |Vn| 1 ' 

is used. 

This is the place to say some words on functionals that include exact exchange. 
As mentioned in Sect. 2.7 this type of functional has been very popular recently 
and improvements of results over GGA-type density functionals for many systems 
and properties have been reported. However, the calculation of the Hartree-Fock 
exchange causes a considerable performance problem in plane wave calculations. 
The Hartree-Fock exchange energy is defined as 604 

where 

Pij (r)=^(r)^(r). (170) 

From this expression the wavefunction force is easily derived and can be calculated 
in Fourier space 

1 dEftpx 



EE V Sfx(G-G')c j (G') . (171) 



The force calculation is best performed in real space, whereas the potential is cal- 
culated in Fourier space. For a system with TVb electronic states and N real space 
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Figure 7. Flow chart for the calculation of the energy and potential of a gradient corrected ex- 
change and correlation functional. 



grid points, a total of 57V, 2 three-dimensional transforms are needed, resulting in 
approximately 25A^ATlogAT operations needed to perform the calculation. This 
has to be compared to the l5Nt,NlogN operations needed for the other Fourier 
transforms of the charge density and the application of the local potential and the 
4N^N operations for the orthogonalization step. In calculations dominated by the 
Fourier transforms an additional factor of at least iVb is needed. If on the other hand 
orthogonalization dominates an increase in computer time by a factor of 5 log N is 
expected. Therefore, at least an order of magnitude more computer time is needed 
for calculations including exact exchange compared to ordinary density functional 
calculations. Consequently, hybrid functional will only be used in exceptional cases 
together with plane waves 262 > 128 . 
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3.4 Total Energy, Gradients, and Stress Tensor 
3.4.1 Total Energy 

Molecular dynamics calculations with interaction potentials derived from density 
functional theory require the evaluation of the total energy and derivatives with 
respect to the parameters of the Lagrangian. In this section formulas are given in 
Fourier space for a periodic system. The total energy can be calculated as a sum of 
kinetic, external (local and non-local pseudopotential) , exchange and correlation, 
and electrostatic energy (to be compared with Eq. (75)) 

£total = E kin + £^ a] + ii'nmilocal + £ xc + E E S ■ (172) 

The individual terms are defined by 



E 



kin 



E^EE^( k )i G + k i 2 i c *( G ' k )i 2 (173) 



2 

G 



^locai -EE A ^cai(G) S/(G)n*(G) (174) 

/ G 

^nTnlocal = £ E /iW £ £ (^iW)* ^^iW (175) 
k i J «,/3e/ 

£ xc = ft£e xc (G)n*(G) (176) 

G 

-Ees = 27rfi £ ! to *^ 2 - + £ ovr i - £ se if. (177) 

G^O 

The overlap between the projectors of the non-local pseudopotential and the Kohn 
Sham orbitals has been introduced in the equation above 

F J ° i (k) = -^£Pi(G)S J (G + k)<(G,k) . (178) 

An alternative expression, using the Kohn-Sham eigenvalues ej(k) can also be used 
£totai = £w k £/i(k)ei(k) 

k i 

-ft£(^xc(G)- £xc (G))n*(G) 

G 

-2uVl £ — J C ^ — + i? ovr i - £ se lf 

G/0 

+AE tot , (179) 

to be compared to Eq. (86). The additional term A£tot m Eq. (179) is needed to 
have an expression for the energy that is quadratic in the variations of the charge 
density, as it is true for Eq. (172). Without the correction term, which is zero for 
the exact charge density, small differences between the computed and the exact 
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density could give rise to large errors in the total energy 129 . The correction energy 
can be calculated from 

G^O ^ ' 

-«E(^"(G) - K°c Ut (G)) (n out (G))*, (180) 

G 

where n m and n out are the input and output charge densities and V™ and V£, ut the 
corresponding exchange and correlation potentials. This term leads to the so-called 
"non-self-consistency correction" of the force, introduced in Eq. (68). 

The use of an appropriate k-point mesh is the most efficient method to calculate 
the total energy of a periodic system. Equivalent, although not as efficient, the 
calculation can be performed using a supcrcell consisting of replications of the 
unit cell and a single integration point for the Brillouin zone. In systems where 
the translational symmetry is broken, e.g. disorder systems, liquids, or thermally 
excited crystals, periodic boundary conditions can still be used if combined with 
a supercell approach. Many systems investigated with the here described method 
fall into these categories, and therfore most calculations using the Car-Parrinello 
molecular dynamics approach are using supercells and a single k-point " integration 
scheme". The only point calculated is the center of the Brillouin zone (r-point; 
k = 0). For the remainder of this chapter, all formulas are given for the r-point 
approximation . 



3.4-2 Wavef unction Gradient 

Analytic derivatives of the total energy with respect to the parameters of the calcu- 
lation are needed for stable molecular dynamics calculations. All derivatives needed 
are easily accessible in the plane wave pseudopotential approach. In the following 
Fourier space formulas are presented 



1 dE tota \ 1 ^ 2 



fi dc*(G) 2 

+ £X(G-G')c j (G') 

G' 

+ J2J2( F uT h ip P p(G)S I (G) , (181) 

where V\ oc is the local potential 

K loc (G) = EA^ cal (G)5 7 (G) + y xc (G)+47r^^ . (182) 

Wavefunction gradients are needed in optimization calculations and in the Car- 
Parrinello molecular dynamics approach. 



58 



3-4-3 Gradient for Nuclear Positions 



The derivative of the total energy with respect to nuclear positions is needed for 
structure optimization and in molecular dynamics, that is 

<9Rz> 9Rj ;S 9R/, S 9R/ ;S 

as the kinetic energy -Ekin and the exchange and correlation energy E xc do not 
depend directly on the atomic positions, the relevant parts are 



dKr.s 



8E PP 

nonlocal 

dKi, s 
dE ES 



-n^iG s A^ cal (G)^(G)n*(G) 



E* E (^)X 
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«E^ ! %P-c(G)^(G) 
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dE OVI \ 



(184) 

<^|(185) 
(186) 



The contribution of the projectors of the non-local pseudopotentials is calculated 
from 



8KZ = -vnV ,G - PliG>s,{G)c '- {G - k) 



(187) 



Finally, the real space part contribution of the Ewald sum is 



dE, 



ovrl 



E'E 




(188) 

The self energy E se \f is independent of the atomic positions and does not contribute 
to the forces. 



3-4-4 Internal Stress Tensor 

For calculations where the supcrccll is changed (e.g. the combination of the Car- 
Parrincllo method with the Parrincllo Rahman approach 201 i 55 ) the electronic in- 
ternal stress tensor is required. The electronic part of the internal stress tensor is 
defined as 440 < 441 (see also Sect. 4.2.3) 

s 

An important identity for the derivation of the stress tensor is 

= n^*)- 1 • (190) 
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The derivatives of the total energy with respect to the components of the cell matrix 
h can be performed on every part of the total energy individually 



ductal dE Mn + dE^ + dE™ iocal | dE xc | dE ES 



dK 



<9h„ 



dh. 



(191) 



Using Eq. (190) extensively the derivatives can be calculated for the case of a plane 
wave basis in Fourier space 202 , 
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(196) 



Finally the derivative of the overlap contribution to the electrostatic energy is 
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(197) 



The local part of the pseudopotential AV^ cal (G) and the nonlocal projector func- 
tions depend on the cell matrix h through the volume, the Bessel transform integral 
and the spherical harmonics function. Their derivatives are lengthy but easy to cal- 
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culate from their definitions Eqs. (140) and (141) 
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(199) 



5.^.5 Non-linear Core Correction 

The success of pseudopotentials in density functional calculations relies on two 
assumptions. The transferability of the core electrons to different environments and 
the linearization of the exchange and correlation energy. The second assumption is 
only valid if the frozen core electrons and the valence state do not overlap. However, 
if there is significant overlap between core and valence densities, the linearization 
will lead to reduced transferability and systematic errors. The most straightforward 
remedy is to include "semi-core states" in addition to the valence shell, i.e. one 
more inner shell (which is from a chemical viewpoint an inert "core level" ) is treated 
cxplicitely. This approach, however, leads to quite hard pseudopotentials which call 
for large plane wave cutoffs. Alternatively, it was proposed to treat the non-linear 
parts of the exchange and correlation energy E xc explicitely 374 . This idea does 
not lead to an increase of the cutoff but ameliorates the above-mentioned problems 
quite a bit. To achieve this, E xc is calculated not from the valence density n(R) 
alone, but from a modified density 

n(R) =n(R) + n coro (R) , (200) 

where n core (R) denotes a density that is equal to the core density of the atomic 
reference state in the region of overlap with the valence density 

n COIC (r) = n COIC (r) if r)r ; (201) 

with the vanishing valence density inside tq. Close to the nuclei a model density 
is chosen in order to reduce the cutoff for the plane wave expansion. Finally, the 
two densities and their derivatives are matched at tq. This procedure leads to a 
modified total energy in Eq. (176), where E xc is replace by 

-Bxc = E xc (n + n corc ) , (202) 

and the corresponding potential is 

V xc = V xc {n + h COTe ) . (203) 



61 



The sum of all modified core densities 

fteoretG) = 5>L. e (G)S/(G) (204) 
I 

depends on the nuclear positions, leading to a new contribution to the forces 



and to the stress tensor 



-fi^;iG 8 V^(G)n^ ore (G)S J (G) , (205) 



^ = EE^c(G)%^5 j( G) . (206) 

The derivative of the core charge with respect to the cell matrix can be performed in 
analogy to the formula given for the local potential. The method of the non-linear 
core correction dramatically improves results on systems with alkali and transition 
metal atoms. For practical applications, one should keep in mind that the non- 
linear core correction should only be applied together with pseudopotentials that 
were generated using the same energy expression. 

3.5 Energy and Force Calculations in Practice 

In Sect. 3.4 formulas for the total energy and forces were given in their Fourier 
space representation. Many terms are in fact calculated most easily in this form, 
but some terms would require double sums over plane waves. In particular, the 
calculation of the charge density and the wavefunction gradient originating from 
the local potential 

XXCG-G'MG') . (207) 

G' 

The expression in Eq. (207) is a convolution and can be calculated efficiently by a 
series of Fourier transforms. The flow charts of this calculations are presented in 
Fig. 8. Both of these modules contain a Fourier transform of the wavefunctions from 
G space to the real space grid. In addition, the calculation of the wavefunction 
forces requires a back transform of the product of the local potential with the 
wavefunctions, performed on the real space grid, to Fourier space. This leads to 
a number of Fourier transforms that is three times the number of states in the 
system. If enough memory is available on the computer the second transform of 
the wavefunctions to the grid can be avoided if the wavefunctions are stored in real 
space during the computation of the density. These modules arc further used in 
the flow chart of the calculation of the local potential in Fig. 9. Additional Fourier 
transforms are needed in this part of the calculation. However, the number of 
transforms does not scale with the number of electrons in the system. Additional 
transforms might be hidden in the module to calculate the exchange and correlation 
potential (see also Fig. 7) and the Poisson solver in cases when the Hockney method 
is used (see Fig. 6). 

The calculation of the total energy, together with the local potential is shown 
in Fig. 10. The overlap between the projectors of the nonlocal pseudopotential and 
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Figure 8. Flow chart for the calculation of the charge density (on the left) and the force on the 
wavefunction from the local potential (on the right). The charge density calculation requires N^, 
(number of states) three dimensional Fourier transforms. For the application of the local potential 
two Fourier transforms per state are needed. If enough memory is available the first transform can 
be avoided if the wavefunction on the real space grid are stored during the density calculation. 
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Figure 9. Flow chart for the calculation of the local potential from the Kohn-Sham orbitals. 
This module calculates also the charge density in real and Fourier space and the exchange and 
correlation energy, Hartrec energy, and local pseudopotcntial energy. 



the wavefunctions calculated in this part will be reused in the calculation of the 
forces on the wavefunctions. There are three initialization steps marked in Fig. 9. 
Step one has only to be performed at the beginning of the calculation, as the 
quantities g and E sc k are constants. The quantities calculated in step two depend 
on the absolute value of the reciprocal space vectors. They have to be recalculated 
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Figure 10. Flow chart for the calculation of the local potential and total energy. Initialization 
steps are marked with numbers. Step 2 has to be repeated whenever the size of the unit cell 
changes. Step 3 has to be repeated whenever nuclear positions have changed. 



whenever the box matrix h changes. Finally, the variables in step three depend 
on the atomic positions and have to be calculated after each change of the nuclear 
positions. The flow charts of the calculation of the forces for the wave functions and 
the nuclei are shown in Figs. 11 and 12. 



3.6 Optimizing the Kohn-Sham Orbitals 

Advances in the application of plane wave based electronic structure methods are 
closely related to improved methods for the solution of the Kohn-Sham equations. 
There are now two different but equally successful approaches available. Fix-point 
methods based on the diagonalization of the Kohn-Sham matrix follow the more 
traditionally ways that go back to the roots of basis set methods in quantum chem- 
istry. Direct nonlinear optimization approaches subject to a constraint were initi- 
ated by the success of the Car-Parrincllo method. The following sections review 
some of these methods, focusing on the special problems related to the plane wave 
basis. 
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Figure 11. Flow chart for the calculation of the forces on the wavefunctions. Notice that the 
calculation of the overlap terms Ff i is done outside the loop over wavefunctions. Besides the 
wavefunctions and the local potential, the structure factors and the projectors of the nonlocal 
pscudopotential are input parameters to this module. 
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Figure 12. Flow chart for the calculation of the forces on the nuclei. 



3.6.1 Initial Guess 

The initial guess of the Kohn Sham orbitals is the first step to a successful calcu- 
lation. One would like to introduce as much knowledge as possible into the first 
step of the calculation, but at the same time the procedure should be general and 
robust. One should also take care not to introduce artifical symmetries that may 
be preserved during the optimization and lead to false results. The most general 
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initialization might be, choosing the wavefunction coefficients from a random dis- 
tribution. It makes sense to weight the random numbers by a function reflecting 
the relative importance of different basis functions. A good choice is a Gaussian 
distribution in G 2 . This initialization scheme avoids symmetry problems but leads 
to energies far off the final results and especially highly tuned optimization methods 
might have problems. 

A more educated guess is to use a superposition of atomic densities and then 
diagonalize the Kohn-Sham matrix in an appropriate basis. This basis can be the 
full plane wave basis or just a part of it, or any other reasonable choice. The 
most natural choice of atomic densities and basis sets for a plane wave calculation 
are the pseudo atomic density and the pseudo atomic wavefunction of the atomic 
reference state used in the generation of the pseudopotcntial. In the CPMD code this 
is one possibility, but often the data needed are not available. For this case the 
default option is to generate a minimal basis out of Slater functions (see Eq. (98) in 
Sect. 2.8) and combine them with the help of atomic occupation numbers (gathered 
using the Aufbau principle) to an atomic density. From the superposition of these 
densities a Kohn-Sham potential is constructed. The Slater orbitals are expanded 
in plane waves and using the same routines as in the standard code the Kohn-Sham 
and overlap matrices arc calculated in this basis. The general eigenvalue problem is 
solved and the cigenfunctions can easily be expressed in the plane wave basis that 
arc in turn used as the initial wavefunctions to the optimization routines. Similarly, 
a given plane wave representation of the total wavefunction can be projected onto an 
auxiliary set of atom-centered functions. This opens up the possibility to perform 
population and bond-order analyses (following for instance the schemes of Mullikcn 
or Mayer) in plane wave-pseudopotential calculations 537 . 



3.6.2 Preconditioning 

Optimizations in many dimensions are often hampered by the appearance of differ- 
ent length scales. The introduction of a metric that brings all degrees of freedom 
onto the same length scale can improve convergence considerably. The applica- 
tion of such a metric is called " preconditioning" and is used in many optimization 
problems. If the variables in the optimization are decoupled the preconditioning 
matrix is diagonal and becomes computationally tractable even for very large sys- 
tems. Fortunately, this is to a large degree the case for a plane wave basis set. For 
large G vectors the Kohn-Sham matrix is dominated by the kinetic energy which 
is diagonal in the plane wave representation. Based on this observation efficient 
preconditioning schemes have been proposed 616 > 61 °j 308 > 344 . The preconditioner im- 
plemented in the CPMD code is based on the diagonal of the Kohn-Sham matrix 
Hg,c, which is given by 

Kg,g' = Hg,g fc.G 1 if|G|>G c (om\ 
K G ,G' = H Gc , Gc 5 G , G , if|G|<G c ' [Zm> 

where G c is a free parameter that can be adjusted to accelerate convergence. How- 
ever, it is found that the actual choice is not very critical and for practical purposes 
it is convenient not to fix G c , but to use an universal constant of 0.5 Hartree for 
-Hg c ,Gc that in turn defines G c for each system. 
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3.6.3 Direct Methods 



The success of the Car-Parrinello approach started the interest in other methods 
for a direct minimization of the Kohn-Sham energy functional. These methods 
optimize the total energy using the gradient derived from the Lagrange function 

C = E KS ({$,}) - ^ «*il*i> - (209) 

ij 

— = ^-^(^1^1^)^. (210) 

3 

Optimization methods differ in the way orbitals are updated. A steepest descent 
based scheme 

c i(G) <— Cj(G) + oKq| g ^(G) (211) 

can be combined with the preconditioner and a line search option to find the optimal 
step size a. Nearly optimal a's can be found with an interpolation based on a 
quadratic polynomial. In Eq. (211) ipi{G) denote the Fourier components of the 
wavefunction gradient. 

Improved convergence can be achieved by replacing the steepest descent step with 
a search direction based on conjugate gradients 594 > 232 > 616 > 23 > 499 

Ci(G) ^ Ci {G) + ahi{G) . (212) 

The conjugate directions are calculated from 
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# ) (G) = K G | G # ) (G) (214) 

< n) = EM n+1) (G)\ 9 ^(G)) 
(^(G)|^(G)) 

A very efficient implementation of this method 616 is based on a band by band 
optimization. A detailed description of this method can also be found in Ref. 472 . 

The direct inversion in the iterative subspace (DIIS) method 495 > 144 > 308 is a 
very successful extrapolation method that can be used in any kind of optimization 
problems. In quantum chemistry the DIIS scheme has been applied to wavefunc- 
tion optimizations, geometry optimizations and in post-Hartree-Fock applications. 
DIIS uses the information of n previous steps. Together with the position vectors 
c-^(G) an estimate of the error vector e^(G) for each previous step k is stored. 
The best approximation to the final solution within the subspace spanned by the 
n stored vectors is obtained in a least square sense by writing 

n 



(G) = J2dkC ( ?\G) , (216) 



fe=i 
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where the dk are subject to the restriction 

n 

]T4 = l 

k=i 

and the estimated error becomes 
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The expansion coefficients dk are calculated from a system of linear equations 
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b kl = J2(eUG)\el(G)) 



(220) 



The error vectors are not known, but can be approximated within a quadratic model 

(221) 



(G) = -K G * #>(G) 



In the same approximation, assuming K to be a constant, the new trial vectors are 
estimated to be 

d(G) = of +1) (G) + K G | G ^" +1) (G) , (222) 
where the first derivative of the energy density functional is estimated to be 

n 

# +1) (G) = ^d fe # ) (G) . 
fe=i 



(223) 



The methods described in this section produce new trail vectors that arc not or- 
thogonal. Therefore an orthogonalization step has to be added before the new 
charge density is calculated 



Cj (G)^5>(G)X, 



(224) 



There are different choices for the rotation matrix X that lead to orthogonal or- 
bitals. Two of the computationally convenient choices are the Lowdin orthogonal- 
ization 



1/2 



and a matrix form of the Gram-Schmidt procedure 



(225) 
(226) 
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where S is the overlap matrix and G is its Cholesky decomposition 

S = GG T . (227) 

Recently new methods that avoid the orthogonalization step have been intro- 
duced. One of them 483 relies on modified functional that can be optimized without 
the orthogonality constraint. These functionals, originally introduced in the con- 
text of linear scaling methods 417 > 452 ; have the property that their minima coincide 
with the original Kohn-Sham energy functional. The methods described above can 
be used to optimize the new functional. 

Another approach 309 is to use a variable transformation from the expansion 
coefficients of the orbitals in plane waves to a set of non-redundant orbital rotation 
angles. This method was introduced in quantum chemistry 618 > 149 > 167 an( j j g usec j 
successfully in many optimization problems that involve a set of orthogonal orbitals. 
A generalization of the orbital rotation scheme allowed the application also for cases 
where the number of basis functions is orders of magnitudes bigger than the number 
of occupied orbitals. However, no advantage is gained over the standard methods, 
as the calculation of the gradient in the transformed variables scales the same as the 
orthogonalization step. In addition, there is no simple and efficient preconditioncr 
available for the orbital rotation coordinates. 



3.6.4 Fix- Point Methods 

Originally all methods to find solutions to the Kohn Sham equations were using 
matrix diagonalization methods. It became quickly clear that direct schemes can 
only be used for very small systems. The storage requirements of the Kohn-Sham 
matrix in the plane wave basis and the scaling proportional to the cube of the basis 
set size lead to unsurmountable problems. Iterative diagonalization schemes can be 
adapted to the special needs of a plane wave basis and when combined with a proper 
preconditioner lead to algorithms that are comparable to the direct methods, both 
in memory requirements and over all scaling properties. Iterative diagonalization 
schemes are abundant. Methods based on the Lanczos algorithm 357 . 151 > 489 can 
be used as well as conjugate gradient techniques 616 ' 97 . Very good results have 
been achieved by the combination of the DIIS method with the minimization of 
the norm of the residual vector 698 ' 344 . The diagonalization methods have to be 
combined with an optimization method for the charge density. Methods based on 
mixing 153 > 4 ; quasi-Newton algorithms 92 > 77 > 319 ) and DIIS 495 > 344 > 345 are successfully 
used. Also these methods use a preconditioning scheme. It was shown that the op- 
timal preconditioning for charge density mixing is connected to the charge dielectric 
matrix 153 > 4 > 299 > 658 > 48 . F or a plane wave basis, the charge dielectric matrix can be 
approximated by expressions very close to the ones used for the preconditioning in 
the direct optimization methods. 

Fix-point methods have a slightly larger prefactor than most of the direct meth- 
ods. Their advantage lies in the robustness and capability of treating systems with 
no or small band gaps. 
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3. 7 Molecular Dynamics 



Numerical methods to integrate the equations of motion are an important part 
of every molecular dynamics program. Therefore an extended literature exists on 
integration techniques (see Ref. 217 and references in there). All considerations 
valid for the integration of equations of motion with classical potentials also apply 
for ab initio molecular dynamics if the Born-Oppenheimer dynamics approach is 
used. These basic techniques will not be discussed here. 

A good initial guess for the Kohn-Sham optimization procedure is a crucial 
ingredient for good performance of the Born-Oppenheimer dynamics approach. An 
extrapolation scheme was devised 24 that makes use of the optimized wavefunctions 
from previous time steps. This procedure has a strong connection to the basic idea 
of the Car-Parrinello method, but is not essential to the method. 

The remainder of this section discusses the integration of the Car-Parrinello 
equations in their simplest form and explains the solution to the constraints equa- 
tion for general geometric constraints. Finally, a special form of the equations of 
motion will be used for optimization purposes. 

3.7.1 Car-Parrinello Equations 

The Car-Parrinello Lagrangian and its derived equations of motions were intro- 
duced in Sect. 2.4. Here Eqs. (41), (44), and (45) are specialized to the case of 
a plane wave basis within Kohn Sham density functional theory. Specifically the 
functions fa are replaced by the expansion coefficients Cj(G) and the orthonormal- 
ity constraint only depends on the wavefunctions, not the nuclear positions. The 
equations of motion for the Car-Parrinello method are derived from this specific 
extended Lagrangian 



where li is the electron mass, and Mi are the masses of the nuclei. Because of 
the expansion of the Kohn Sham orbitals in plane waves, the orthonormality con- 
straint does not depend on the nuclear positions. For basis sets that depend on 
the atomic positions (e.g. atomic orbital basis sets) or methods that introduce an 
atomic position dependent metric (ultra-soft pseudopotentials 661 > 351 ; PAW 143 > 347 ; 
the integration methods have to be adapted (see also Sect. 2.5). Solutions that in- 
clude these cases can be found in the literature 280 > 351 . 143 < 310 . The Euler-Lagrangc 
equations derived from Eq.( 228) are 
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The two sets of equations are coupled through the Kohn-Sham energy functional 
and special care has to be taken for the integration because of the orthonormality 
constraint. 

The integrator used in the CPMD code is based on the velocity Verlet / RATTLE 
algorithm 603 ! 638 > 15 . The velocity Verlet algorithm requires more operations and 
more storage than the Verlet algorithm 664 . However, it is much easier to incorpo- 
rate temperature control via velocity scaling into the velocity Verlet algorithm. In 
addition, velocity Verlet allows to change the time step trivially and is conceptually 
easier to handle 638 > 391 . It is defined by the following equations 



R I (t + St) = R I (t) + ^-F I (t) (231) 

R/(t + St) = Rj(t) + St Ri{t + St) 
t I (t + 5t) = c I (t) + -^f i (t) 
Ci(t + St) = a(t) + St 5i(t + St) 

a(t + st) = c t {t + st) + ^2 x y c j(*) 

3 

calculate F/(i+£i) 
calculate f, (t + St) 

Rj(t + St) = Ri(t + St) + F/(£ + St) 
c'i{t + St) = ii{t + St) + -^fi(i + St) 

C, (t + St) = C , l (t+St)+J2 Y *3 C 3 (* + <**) . 



where R/ (t) and c j (t) are the atomic positions of particle / and the Kohn-Sham 
orbital i at time t respectively. Here, Fj arc the forces on atom /, and f$ are the 
forces on Kohn-Sham orbital i. The matrices X and Y are directly related to the 
Lagrange multipliers by 



X ^ - 2^ 



V — 



St_ 



Ay, 



(232) 
(233) 



Notice that in the rattle algorithm the Lagrange multipliers to enforce the or- 
thonormality for the positions A p and velocities A v are treated as independent 
variables. Denoting with C the matrix of wavefunction coefficients Cj(G), the or- 
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thonormality constraint can be written as 

&(t + 6t)C(t + 6t) -1 = (234) 

C + Xcj 1 fc + Xcj -1 = (235) 

C t C + XC t C + C t CX t +XX f -I = (236) 
XX f + XB + B t X t = I A , (237) 

where the new matrices Ay = c- (t + St)cj(t + 5t) and By = c\ (t)cj(t + 6t) have 
been introduced in Eq. (237). The unit matrix is denoted by the symbol I. By 
noting that A = I + 0(St 2 ) and B = I + O(St), Eq. (237) can be solved iteratively 
using 

1 



x (n+i) = ± 
2 



I A + X (n) (I - B) 
+ (I-B)XW - (x^) 2 



and starting from the initial guess 

X (0) 



1 



(I-A) . 



(238) 



(239) 



In Eq. (238) it has been made use of the fact that the matrices X and B are real 
and symmetric, which follows directly from their definitions. Eq. (238) can usually 
be iterated to a tolerance of 10 -6 within a few iterations. 

The rotation matrix Y is calculated from the orthogonality condition on the 
orbital velocities 

c\ (t + 5t)cj (t + St) + c\(t + 5t)c 3 (t + 6t)= 0. (240) 
Applying Eq. (240) to the trial states C' + YC yields a simple equation for Y 



(241) 



where Q y = c\{t + 5t)c'\{t+5t). The fact that Y can be obtained without iteration 
means that the velocity constraint condition Eq. (240) is satisfied exactly at each 
time step. 



3.7.2 Advanced Techniques 

One advantage of the velocity Verlet integrator is that it can be easily combined 
with multiple time scale algorithms 636 > 639 and still results in reversible dynamics. 
The most successful implementation of a multiple time scale scheme in connec- 
tion with the plane wavc-pseudopotcntial method is the harmonic reference system 
idea 471 > 639 . The high frequency motion of the plane waves with large kinetic energy 
is used as a reference system for the integration. The dynamics of this reference 
system is harmonic and can be integrated analytically. In addition, this can be com- 
bined with the basic notion of a prcconditioner already introduced in the section 
on optimizations. The electronic mass used in the Car-Parrinello scheme is a ficti- 
tious construct (see Sect. 2.4, Eq. (45)) and it is allowed to generalize the idea by 
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introducing different masses for different "classical" degrees of freedom 473 ,6io,639^ 
In agreement with the preconditioner introduced in the optimization section, the 
new plane wave dependent masses are 

M(G) = {(n/a) (±G 2 +V(G,G)) H(G^ G)> a ' (242) 

where H and V are the matrix elements of the Kohn-Sham matrix and the poten- 
tial respectively. The reference electron mass is fi and the parameter a has been 
introduced before in Eq. (208) as Hq c ,g c - With the preconditioned masses and 
the harmonic reference system, the equations of motion of the system are 

/i(G)c j (G) = -A(G)c i (G) + *$j(G) + 5^Ayc J -(G) . (243) 

3 

where 5$j(G) is the force on orbital i minus — A(G). From Eq. (243) it is easy 
to see that the frequencies u(G) = ^/A(G)//u(G) are independent of G and that 
there is only one harmonic frequency equal to wa/fl. The revised formulas for 
the integration of the equations of motion for the velocity Verlet algorithm can be 
found in the literature 639 . 

The implications of the G vector dependent masses can be seen by revisiting 
the formulas for the characteristic frequencies of the electronic system Eqs. (52), 
(53), and (54). The masses fi arc chosen such that all frequencies u>ij are approxi- 
mately the same, thus optimizing both, adiabaticity and maximal time step. The 
disadvantage of this method is that the average electron mass seen by the nuclei is 
drastically enhanced, leading to renormalization corrections on the masses Mi 
that are significantly higher than in the standard approach and not as simple to 
estimate by an analytical expression. 

3.7.3 Geometrical Constraints 

Geometrical constraints are used in classical simulations to freeze fast degrees of 
freedom in order to allow for larger time steps. Mainly distance constraints are 
used for instance to fix intramolecular covalent bonds. These type of applications 
of constraints is of lesser importance in ab initio molecular dynamics. However, in 
the simulation of rare events such as many reactions, constraints play an important 
role together with the method of thermodynamic integration 217 . The "blue-moon" 
ensemble method 115 » 589 enables one to compute the potential of mean force. This 
potential can be obtained directly from the average force of constraint and a geo- 
metric correction term during a molecular dynamics simulation as follows: 

Hb)-HZi) = / 5 ^'(-^f) ' (244) 

where T is the free energy and £(r) a one-dimensional reaction coordinate, 7i the 
Hamiltonian of the system and (• • • }|° nd - the conditioned average in the constraint 
ensemble 589 . By way of the blue moon ensemble, the statistical average is replaced 
by a time average over a constrained trajectory with the reaction coordinate fixed 
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at special values, £(R) = £' , and £(R, R) = 0. The quantity to evaluate is the 
mean force 



dT (Z-V2 hA + fc B TG]) e , 



(245) 



(z- 1/2 ) e 

where A is the Lagrange multiplier of the constraint, 



and 

Z 2 ^ MjMj 9Rj 9R/9R,/ 5Rj ' 1 j 

where (• • • is the unconditioned average, as directly obtained from a constrained 
molecular dynamics run with £(R) = £' an d 

m)-HZi) = J <%'— (248) 

finally defines the free energy difference. For the special case of a simple distance 
constraint £(R) = |R/ — Rj| the parameter Z is a constant and G = 0. 

The rattle algorithm, allows for the calculation of the Lagrange multiplier of 
arbitrary constraints on geometrical variables within the velocity Verlet integrator. 
The following algorithm is implemented in the CPMD code. The constraints are 
defined by 

<7«({R / (i)}) = , (249) 
and the velocity Verlet algorithm can be performed with the following steps. 

R 7 = R/(f) + 6tRi 
Rtf + St) =R J + A_ gp (t) 
calculate Fi(t + St) 

R'^Ri + ^-F^t + St) 

Ri(t + St) = R'/ + ^gv(t + St) , 

where the constraint forces are defined by 

da^({Hj(t)}) 



8p(*) = - 2^ A p ( 25 °) 

SvW = -2^ A V QJ^ • ( 251 ) 
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The Lagrange multiplier have to be determined to ensure that the constraint on 
the positions and velocities are exactly fulfilled at the end of the time step. For the 
position, the constraint condition is 

a w ({R/(t + ^)}) = . (252) 

Eq. (252) is in general a system of nonlinear equations in the Lagrange multipliers 
Ap. These equations can be solved using a generalized Newton algorithm 491 that 
can be combined with a convergence acceleration scheme based on the direct inver- 
sion in the iterative subspace method 495 ' 144 . The error vectors for a given set of 
Lagrange multipliers A are calculated from 

e J (A) = -^J J -- 1 (A) ( r«(A) . (253) 

o 

The Jacobian J is defined by 

1/2 f/WKO) , (256) 



= -E 



2M 



where ff(A) = XL \ l da^ /dRi. Typically only a few iterations are needed to 
converge the Lagrange multipliers to an accuracy of 1 x 1CP 8 . 

The constraint condition for the velocities can be cast into a system of linear 
equations. Again, as in the case of the orthonormality constraints in the Car- 
Parrinello method, the Lagrange multiplier for the velocity update can be calculated 
exactly without making use of an iterative scheme. Defining the derivative matrix 



the velocity constraints are 



8(7^ 

A„ = w , (257) 



& {i \t + 5t)=Q (258) 



E 



R/ = (259) 



- E (E ^ A " A ^ J A J = E ■ (260) 

The only information needed to implement a new type of constraint are the formulas 
for the functional value and its derivative with respect to the nuclear coordinates 
involved in the constraint. 
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3. 7.4 Using Car-Parrinello Dynamics for Optimizations 



By adding a friction term, Car-Parrincllo molecular dynamics can be turned into 
a damped second order dynamics scheme (see also Sect. 2.4.6). 

The friction can be applied both to the nuclear degrees of freedom and the 
electronic coordinates. The resulting dynamics equation are a powerful method to 
simultaneously optimize the atomic structure and the Kohn-Sham orbitals 472 > 610 . 
Harmonic reference system integration and plane wave dependent electron masses, 
introduced above, are especially helpful in this context, as the derived dynamics 
does not have a direct physical relevance. 

Introducing a friction force proportional to the constants 7„ and 7 e the equations 
of motion can readily be integrated using the velocity Vcrlet algorithm. The friction 
terms translate into a simple rescaling of the velocities at the beginning and end of 
the time step according to 

R/(t) = 7nR/(t) 

Ci(t) = 7 c c t (i) 
VELOCITY VERLET UPDATE 

K T {t+6t) = 7n R 7 (t + (5t) 
Ci(t+6t) = j e Ci{t + St) . 

It was shown 472 ' 610 that this scheme leads to optimizations that are competitive 
with other methods described in Sect. 3.6 



3.8 Data Structures and Computational Kernels 

In the practical implementation of the method, mathematical symbols have to 
be translated into data structures of the computer language. Then mathematical 
formulas are set into computer code using the data structures. The layout of the 
data structures should be such that optimal performance for the algorithms can be 
achieved. The CPMD code is written in FORTRAN77 and in the following sections the 
most important data structures and computational kernels will be given in pseudo 
code form. The following variables are used to denote quantities that measure 



system size. 

N at number of atoms 

-/Vp number of projectors 

TVb number of electronic bands or states 

TVpw number of plane-waves 

TVd number of plane-waves for densities and potentials 
N x , N y , N z number of grid points in x, y, and z direction 



N = N K N y N z total number of grid points 

In Table 3 the relative size of this variables are given for two systems. The example 
for a silicon crystal assumes an energy cutoff of 13 Rydberg and s non-locality for 
the pscudopotcntial. In the example of a water system the numbers are given per 
molecule. The cutoff used was 70 Rydberg and the oxygen pseudopotential has a s 
nonlocal part, the hydrogen pseudopotential is local. 
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Table 3. Relative size of characteristic variables in a plane wave calculation. See text for details. 





silicon 


water 




1 


3 


^P 


1 


1 


N h 


2 


4 


AW 


53 


1000 


A^D 


429 


8000 


N 


1728 


31250 



3.8.1 CPMD Program: Data Structures 

Important quantities in the pscudopotcntial plane wave method depend cither not 
at all, linearly, or quadratically on the system size. Examples for the first kind of 
data are the unit cell matrix h and the cutoff E cut . Variables with a size that grows 
linearly with the system are 

T(3,N at ) nuclear positions 

v(3,N at ) nuclear velocities 

f ( 3 , N at ) nuclear forces 

g ( 3 , Npw ) plane-wave indices 

ipg(3, Npw) mapping of G-vectors (positive part) 

img(3, Npw) mapping of G-vectors (negative part) 

r~h.og(Npw) densities (n, n c , n to t) in Fourier-space 

vpot (iVpw) potentials (Vl oc , V xc , Vh) in Fourier-space 

n(N x ,N y ,N Z ) densities (n, n c , n to t) in real-space 

v(N x ,N y ,N Z ) potentials (Vioc, V xc , Vh) in real-space 

vps(Nrj) local pseudopotential 

rpc(A£>) core charges 

pro (Npw) projectors of non-local pseudopotential. 

The pseudopotential related quantities vps, rpc, and pro are one-dimensional in 
system size but also depend on the number of different atomic species. In the 
following it is assumed that this is one. It is easy to generalize the pseudo codes 
given to more than one atomic species. For real quantities that depend on G— 
vectors only half of the values have to be stored. The other half can be recomputed 
when needed by using the symmetry relation 



A(G) = A*(-G) . (261) 

This saves a factor of two in memory. In addition G vectors are stored in a linear 
array, instead of a three-dimensional structure. This allows to store only non-zero 
variables. Because there is a spherical cutoff, another reduction of a factor of two 
is achieved for the memory. For the Fourier transforms the variables have to be 
prepared in a three-dimensional array. The mapping of the linear array to this 
structure is provided by the information stored in the arrays ipg and img. 
Most of the memory is needed for the storage of quantities that grow quadratically 
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with system size. 



eigr ( Np , N at ) structure factors 

fnl(N p ,Nb) overlap of projectors and bands 

df nl (N p ,N b ,3) derivative of f nl 

smat (Nb, Nb) overlap matrices between bands 

cr {Npw , Nb) bands in Fourier space 

cv(Npw ,Nb) velocity of bands in Fourier space 

cf (Npw ,Nb) forces of bands in Fourier space 



In order to save memory it is possible to store the structure factors only for the 
G vectors of the wave function basis or even not to store them at all. However, 
this requires that the missing structure factors are recomputed whenever needed. 
The structure factors eigr and the wavefunction related quantities cr , cv, cf are 
complex numbers. Other quantities, like the local pseudopotential vps, the core 
charges rpc, and the projectors pro can be stored as real numbers if the factor 
(— i) 1 is excluded. 

3.8.2 CPMD Program: Computational Kernels 

Most of the calculations in a plane wave code are done in only a few kernel routines. 
These routines are given in this section using a pseudo code language. Where 
possible an implementation using basic linear algebra (blas) routines is given. The 
first kernel is the calculation of the structure factors. The exponential function of 
the structure factor separates in three parts along the directions s x , s y , s z . 

MODULE StructureFactor 
FOR i=l:iV at 

s(l:3) = 2 * PI * MATMUL [html (1:3, 1:3) ,r(l:3,i)] 
dp(l:3) = CMPLX[C0S[s(l:3)] , SIN [s (1:3)]] 
dm(l:3) = C0NJG[dp(l :3)] 
e(0,l:3,i) = 1 
FOR k=l: ffmax 

e(k,l:3,i) = e(k-l,l:3,i) * dp 

e(-k,l:3,i) = e (-k+1 , 1 : 3 , i) * dm 

END 

FOR j=0:N B 

eigr(j,i) = e(g(l,j) ,l,i) * e (g(2 , j ) , 2 , i) * e(g(3, j) ,3,i) 

END 



In the module above html is the matrix (h*) -1 . One of the most important calcu- 
lation is the inner product of two vectors in Fourier space. This kernel appears for 
example in the calculation of energies 



Making use of the fact that both functions are real the sum can be restricted to half 



END 




(262) 



G 
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of the G vectors, and only real operations have to be performed. Approximately a 
factor of three in operations can be save this way. Special care has to be taken for 
the zero G vector. It is assumed that this plane wave component is stored in the 
first position of the arrays. 

MODULE DotProduct 
e = A(l) * B(l) 
FOR i=2:iV D 

ar = REAL (A (i)) 

ai = IMAG(A(i) ) 

br = REAL(B(i)) 

bi = IMAG(B(i)) 

e = e + 2 * (ar * br + ai * bi) 



This loop structure is available in the blas library, optimized on most computer 
architectures. To use the BLAS routines for real variables, complex numbers have 
to be stored as two real numbers in contiguous memory locations. 

e = A(l) * B(l) + 2 * sdot(2 * iV D - 2,A(2) , 1 ,B(2) , 1) 

The calculation of overlap matrices between sets of vectors in real space is a im- 
portant task in the orthogonalization step 



It can be executed by using matrix multiply routines from the BLAS library. The 
special case of the zero G vector is handled by a routine that performs a rank 1 
update of the final matrix. 

MODULE Overlap 

CALL SGEMM ( ' T ' , 'N' ,N h ,N h ,2*N PW ,2,t 

& ca(l,l) ,2*AW,cb(l,l) ,2*iVp W ,0,smat,7Vb) 
CALL SDER(7V b ,^Vb,-l,ca(l,l) ,2*iVp W ,cb(l,l) ,2*iVp W ,smat,7Vb) 

For a symmetric overlap additional time can be saved by using the symmetric 
matrix multiply routine. The overlap routines scale like N^Npyj- It is therefore 
very important to have an implementation of these parts that performs close to 
peek performance. 

MODULE SymmetricOverlap 

CALL SSYRKCU' , 'T' , N b , 2*iV PW , 2 , ca(l , 1) , 2*N PW , , smat , N h ) 
CALL SDER(7V b ,^Vb,-l,ca(l,l) ,2*iVp W ,cb(l,l) ,2*iVp W ,smat,7Vb) 

Another operation that scales as the overlap matrix calculations is the rotation of 



END 




(263) 



G 
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a set of wavefunctions in Fourier space 

B i (G) = J2MG)S ji . (264) 

3 

Again this kernel can be executed by using the optimized blas matrix multiply 
routines. 

MODULE Rotation 

CALL SGEMM ( ' N ' , >N' , 2*iV PW , N h , , 1 , ca( 1 , 1) , 2*iV PW ,& 
& smat,A^ b ,0,cb(l,l) ,2*iVp W ) 

The overlap calculation of the projectors of the nonlinear pseudopotential with 
the wavefunctions in Fourier space scales as N p N^Npyj. As the projectors are 
stored as real quantities, the imaginary prefactor and the structure factor have 
to be applied before the inner product can be calculated. The following pseudo 
code calculates M projectors at a time, making use of the special structure of the 
prefactor. This allows again to do all calculations with real quantities. The code 
assumes that the total number of projectors is a multiple of M . A generalization 
of the code to other cases is straightforward. By using batches of projectors the 
overlap can be calculated using matrix multiplies. The variable lp(i) holds the 
angular momentum of projector i. 

MODULE FNL 
FOR i=l:iV p ,M 

IF (M0D(lp(i) ,2) == 0) THEN 
FOR j=0:M-l 

pf = -l**(lp(i+j)/2) 
FOR k=l:iVpw 

t = pro(k) * pf 
er = REAL[eigr(k,iat(i+j))] 
ei = IMAG[eigr(k,iat(i+j))] 
scr(k,j) = CMPLX[t * er,t * ei] 

END 

END 
ELSE 

FOR j=0:M-l 

pf = -l**(lp(i+j)/2+l) 
FOR k=l:iV PW 

t = pro(k) * pf 
er = REAL[eigr(k,iat(i+j))] 
ei = IMAG[eigr(k,iat(i+j))] 
scr(k,j) = CMPLX[-t * ei,t * er] 

END 

END 
END IF 

scr(l,0:M-l) = scr(l,0:M-l)/2 
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CALL SGEMM ( ' T ' , 'N' ,M, N h ,2*N PW ,2 ,k 

& scr(l,0) ,2*iVp W ,cr(l,l) ,2*A r Pw.0,fnl(i,l) ,7V p ) 

END 

Fourier transform routines are assumed to work on complex data and return 
also arrays with complex numbers. The transform of data with the density cutoff is 
shown in the next two pseudo code sections. It is assumed that a three dimensional 
fast Fourier transform routine exists. This is in fact the case on most computers 
where optimized scientific libraries are available. The next two pseudo code seg- 
ments show the transform of the charge density from Fourier space to real space 
and back. 

MODULE INVFFT 
scr(l:A^ x ,l:iV y ,l:Ar z ) = 
FOR i=l:N B 

scr(ipg(l,i) ,ipg(2,i) ,ipg(3,i)) = rhog(i) 

scr(img(l,i) ,img(2,i) ,img(3,i)) = CONJG[rhog(i)] 

END 

CALL FFT3D("INV",scr) 

n(l:N x ,l:N y ,l:N z ) = REAL [scr (1 : N x , 1 : N y , 1 : N z )~] 
MODULE FWFFT 

scr(l:iV x ,l:iV y> l:Ar z ) = n(l : JV X , 1 : N y , 1 : N z ) 
CALL FFT3D("FW",scr) 
FOR i=l:iV D 

rhog(i) = scr(ipg(l,i) ,ipg(2,i) ,ipg(3,i)) 

END 

Special kernels are presented for the calculation of the density and the application 
of the local potential. These are the implementation of the flow charts shown in 
Fig. 8. The operation count of these routines is N\>Nlog[N]. In most applications 
these routines take most of the computer time. Only for the biggest applications 
possible on todays computers the cubic scaling of the orthogonalization and the 
nonlocal pseudopotential become dominant. A small prefactor and the optimized 
implementation of the overlap are the reasons for this. 

In the Fourier transforms of the wavefunction two properties are used to speed 
up the calculation. First, because the wavefunctions are real two transforms can be 
done at the same time, and second, the smaller cutoff of the wavefunctions can be 
used to avoid some parts of the transforms. The use of the sparsity in the Fourier 
transforms is not shown in the following modules. In an actual implementation a 
mask will be generated and only transforms allowed by this mask will be executed. 
Under optimal circumstances a gain of almost a factor of two can be achieved. 
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MODULE Density 

rho(l:7V x ,l:iV y ,l:./V z ) = 
FOR i=l:JV b ,2 

scT(l:N K ,l:N y ,l:N z ) = 

FOR j = l:7V PW 

scr(ipg(l,i) ,ipg(2,i) ,ipg(3,i)) = c(j,i) + I * c(j,i+l) 
scr(img(l,i) ,img(2,i) ,img(3,i)) = CONJG[c(j,i) + I * c(j,i+l)] 

END 

CALL FFT3D("INV",scr) 

rho(l:Ar x ,l:JV y> l:Ar z ) = rho(l : N K , 1 : N y , 1 : N z ) + t 

& REAL [scr (l:./V x ,l:iV y ,l:iV z )]**2 + IMAG [scr (1 : N x , 1 : N y , 1 : AT Z )] **2 

END 

MODULE VPSI 
FOR i=l:iV b ,2 

scr(l:N x ,l:N y ,l:N z ) = 

FOR j=l:iVp W 

scr(ipg(l,i) ,ipg(2,i) ,ipg(3,i)) = c(j,i) + I * c(j,i+l) 
scr(img(l,i) ,img(2,i) ,img(3,i)) = CONJG[c(j,i) + I * c(j,i+l)] 

END 

CALL FFT3D("INV",scr) 

scr(l:iV x ,l:Ar y ,l:iV z ) = scr (1 : N K , 1 : N y , 1 : N z ) * t 

& vpot(l:JV x> l:Ar y ,l:iV z ) 
CALL FFT3D("FW",scr) 
FOR j=l:7V PW 

FP = scr(ipg(l,i) ,ipg(2,i) ,ipg(3,i)) & 

& + scr (irngCl , i) ,img(2,i) ,img(3,i)) 

FM = scr(ipg(l,i) ,ipg(2,i) ,ipg(3,i)) & 

& - scr(img(l,i) ,img(2,i) ,img(3,i)) 

fc(j,i) = f(i) * CMPLX [REAL [FP] , IMAG [FM] ] 

fc(j,i+l) = f(i+l) * CMPLX [IMAG [FP] , -REAL [FM] ] 

END 

END 



3.9 Parallel Computing 
3. 9. 1 Introduction 

Ab initio molecular dynamics calculation need large computer resources. Memory 
and CPU time requirement make it necessary to run projects on the biggest com- 
puters available. It is exclusively parallel computers that provide these resources 
today. There are many different types of parallel computers available. Comput- 
ers differ in their memory access system and their communication system. Widely 
different performances are seen for bandwidth and latency. In addition, different 
programming paradigms are supported. In order to have a portable code that can 
be used on most of the current computer architectures, CPMD was programmed us- 
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ing standard communication libraries and making no assumption on the topology 
of the processor network and memory access system. 

Minimizing the communication was the major goal in the implementation of 
the parallel plane wave code in CPMD. Therefore, the algorithms had to be adapted 
to the distributed data model chosen. The most important decisions concern the 
data distribution of the largest arrays in the calculation. These arrays are the ones 
holding information on the wavefunctions. Three distribution strategies can be 
envisaged and were used before 9( M37,687,688,m^ 

First, the data are distributed over the bands 687 . Each processor holds all 
expansion coefficients of an electronic band locally. Several problems arise with 
this choice. The number of bands is usually of the same magnitude as the number 
of processors. This leads to a severe load-balancing problem that can only be 
avoided for certain magic numbers, namely if the number of bands is a multiple 
of the number of CPU's. Furthermore this approach requires to perform three- 
dimensional Fourier transforms locally. The memory requirements for the Fourier 
transform only increase linearly with system size, but their prefactor is very big and 
a distribution of these arrays is desirable. In addition, all parts of the program that 
do not contain loops over the number of bands have to be parallelized using another 
scheme, leading to additional communication and synchronization overhead. 

Second, the data is distributed over the Fourier space components and the 
real space grid is also distributed 90 > 137 > 117 . This scheme allows for a straight 
forward parallclization of all parts of the program that involve loops over the Fourier 
components or the real space grid. Only a few routines are not covered by this 
scheme. The disadvantage is that all three-dimensional Fourier transforms require 
communication. 

Third, it is possible to use a combination of the above two schemes 688 . This 
leads to the most complicated scheme, as only a careful arrangement of algorithms 
avoids the disadvantages of the other schemes while still keeping their advantages. 

Additionally, it is possible to distribute the loop over k-points. As most calcu- 
lation only use a limited number of k-points or even only the T-point, this method 
is of limited use. However, combining the distribution of the k-points with one of 
the other method mentioned above might result in a very efficient approach. 

The CPMD program is parallelized using the distribution in Fourier and real space. 
The data distribution is held fixed during a calculation, i.e. static load balancing 
is used. In all parts of the program where the distribution of the plane waves does 
not apply, an additional parallclization over the number of atoms or bands is used. 
However, the data structures involved are replicated on all processors. 

A special situation exists for the case of path integral calculations (see Sect. 4.4), 
where an inherent parallclization over the Trotter slices is present. The problem is 
"embarrassingly parallel" in this variable and perfect parallelism can be observed on 
all types of computers, even on clusters of workstations or supercomputers ("meta- 
computing" ) . In practice the parallclization over the Trotter slices will be combined 
with one of the schemes mentioned above, allowing for good results even on mas- 
sively parallel machines with several hundred processors. 
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3.9.2 CPMD Program: Data Structures 



In addition to the variables used in the serial version, local copies have to be de- 
fined. These local variables will be indexed by a superscript indicating the processor 
number. The total number of processors is P. Each processor has a certain number 
of plane waves, atoms, electronic bands and real space grid points assigned. 

N% t number of atoms on processor p 

NP number of projectors on processor p 

N£ number of electronic bands or states on processor p 

Np W number of plane- waves on processor p 

iVp number of plane- waves for densities and potentials on processor p 

N^, N y , N z number of grid points in x, y, and z direction on processor p 
N p —NPN y N z total number of grid points on processor p 

The real space grid is only distributed over the x coordinates. This decision is 
related to the performance of the Fourier transform that will be discussed in more 
detail in the following sections. The distribution algorithm for atoms, projectors 
and bands just divides the total number of these quantities in equal junks based on 
their arbitrary numbering. The algorithms that use these parallclization schemes 
do not play a major role in the overall performance of the program (at least for the 
systems accessible with the computers available today) and small imperfections in 
load balancing can be ignored. 

Data structures that are replicated on all processors: 
r(3,N at ) nuclear positions 

v(3,N at ) nuclear velocities 

i(3,N at ) nuclear forces 

inl(N p ,N},) overlap of projectors and bands 
smat (Nb,Ni,) overlap matrices between bands. 

Data structures that are distributed over all processors: 

g ( 3 , Np W ) plane- wave indices 

ipg(3, Np W ) mapping of G-vectors (positive part) 

±mg(3,N PW ) mapping of G-vectors (negative part) 

rhogCiVp^) densities (n, n c , n to t) in Fourier-space 

vpot(iVp W ) potentials (Vi oc , Vx C , Vh) in Fourier-space 

tl(NP ,N y ,N Z ) densities (n, n c , ntot) in real-space 

v(NP ,N y ,N Z ) potentials (Vioc, Vx C , Vh) in real-space 

vpsCiV^,) local pseudopotential 

rpc(iV^) core charges 

pro (N pw ) projectors of non-local pseudopotential 

eigr ( , N at ) structure factors 

df nl (7V p , Nl , 3) derivative of f nl 

cr (Np W , Ni,) bands in Fourier space 

cv(N PW , Nb) velocity of bands in Fourier space 

cf ( N PW , Nb) forces of bands in Fourier space. 

Several different goals should be achieved in the distribution of the plane waves 
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over processors. All processors should hold approximately the same number of 
plane waves. If a plane wave for the wavefunction cutoff is on a certain processor, 
the same plane wave should be on the same processor for the density cutoff. The 
distribution of the plane waves should be such that at the beginning or end of a 
three dimensional Fourier transform no additional communication is needed. To 
achieve all of these goals the following heuristic: algorithm 137 is used. The plane 
waves are ordered into "pencils". Each pencil holds all plane waves with the same 
g y and g z components. The pencils are numbered according to the total number 
of plane waves that are part of it. Pencils are distributed over processors in a 
"round robin" fashion switching directions after each round. This is first done for 
the wavefunction cutoff. For the density cutoff the distribution is carried over, and 
all new pencils are distributed according to the same algorithm. Experience shows 
that this algorithm leads to good results for the load balancing on both levels, the 
total number of plane waves and the total number of pencils. The number of pencils 
on a processor is proportional to the work for the first step in the three-dimensional 
Fourier transform. 

Special care has to be taken for the processor that holds the G = component. 
This component has to be treated individually in the calculation of the overlaps. 
The processor that holds this component will be called pO. 



3.9.3 CPMD Program: Computational Kernels 

There arc three communication routines mostly used in the parallelization of the 
CPMD code. All of them are collective communication routines, meaning that all 
processors are involved. This also implies that synchronization steps are performed 
during the execution of these routines. Occasionally other communication routines 
have to be used (e.g. in the output routines for the collection of data) but they 
do not appear in the basic computational kernels. The three routines are the 
Broadcast, GlobalSum, and MatrixTranspose. In the Broadcast routine data is 
send from one processor (px) to all other processors 

x p <- a; px . (265) 

In the GlobalSum routine a data item is replaced on each processor by the sum over 
this quantity on all processors 

x p ^j2 xP ■ ( 266 ) 

p 

The MatrixTranspose changes the distribution pattern of a matrix, e.g. from row 
distribution to column distribution 

x(p,:)^x(:,p) . (267) 

On a parallel computer with P processors, a typical latency time t^ (time for the 
first data to arrive) and a bandwidth of B, the time spend in the communication 
routines is 

Broadcast log 2 [P] {t L + N/B} 

GlobalSum log 2 [P] {t L + N/B} 

MatrixTranspose Pt L + N/(PB) 
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Table 4. Distribution of plane waves and "pencils" in parallel runs on different numbers of pro- 
cessors. Example for a cubic box with a volume of 6479.0979 bohr 3 and a 70 Rydberg cutoff for 
the wavefunctions. This is the simulation box needed for 32 water molecules at normal pressure. 
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where it is assumed that the amount of data N is constant. The time needed 
in Broadcast and GlobalSum will increase with the logarithm of the number of 
processors involved. The time for the matrix transposition scales for one part 
linearly with the number of processors. Once this part is small, then the latency 
part will be dominant and increase linearly. Besides load balancing problems, the 
communication routines will limit the maximum speedup that can be achieved on 
a parallel computer for a given problem size. Examples will be shown in the last 
part of this section. 

With the distribution of the data structures given, the parallelization of the com- 
putational kernels is in most cases easy. In the StructureFactor and Rotation 
routines the loop over the plane waves iNfo has to be replaced by . The routines 
performing inner products have to be adapted for the G = term and the global 
summation of the final result. 

MODULE DotProduct 
IF (p == P0) THEN 

ab = A(l) * B(l) + 2 * sdot(2 * (iVg - 1) , A(2) , 1 ,B(2) , 1) 
ELSE 
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ab = 2 * sdot(2 * , A(l) , 1 ,B(1) , 1) 
END IF 

CALL GlobalSum[ab] 



MODULE Overlap 

CALL SGEMM ( ' T ' , 'N' ,N h , N h ,2*iV£ w ,2, & 

& ca(l,l) ,2*iV^ w> cb(l,l) ,2*N]p W ,0,sma.t,N b ) 
IF (p == PO) CALL SDERUV b ,iV b ,-l,ca(l,l),2*iV£ w> & 

& cb(l,l) ,2*A r ^ w ,smat,iVb) 
CALL GlobalSum[smat] 



Similarly, the overlap part of the FNL routine has to be changed and the loops 
restricted to the local number of plane waves. 



MODULE FNL 
FOR i=l:iV p ,M 

IF (M0D(lp(i) ,2) == 0) THEN 
FOR j=0:M-l 

pf = -l**(lp(i+j)/2) 
FOR k=l:7V£ w 

t = pro(k) * pf 
er = REAL[eigr(k,iat(i+j))] 
ei = IMAG[eigr(k,iat(i+j))] 
scr(k,j) = CMPLX[t * er,t * ei] 

END 

END 
ELSE 

FOR j=0:M-l 

pf = -l**(lp(i+j)/2+l) 
FOR k=l:7V£ w 

t = pro(k) * pf 
er = REAL[eigr(k,iat(i+j))] 
ei = IMAG[eigr(k,iat(i+j))] 
scr(k,j) = CMPLX[-t * ei,t * er] 

END 

END 
END IF 

IF (p == PO) scr(l,0:M-l) = scr(l,0:M-l)/2 
CALL SGEMM ( ' T ' , >N' ,M, N h , 2*N? W , 2 ,& 

k scr(l,0),2*iVP w ,cr(l,l),2*^ w ,0,fnl(i,l),7V p ) 

END 

CALL GlobalSum[fnl] 



The routines that need the most changes are the once that include Fourier trans- 
forms. Due to the complicated break up of the plane waves a new mapping has to 
be introduced. The map mapxy ensures that all pencils occupy contiguous memory 
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locations on each processor. 



MODULE INVFFT 

scrl(l:iV x> l:^ p D oncil ) = 
FOR i=l:iVg 

scrl(ipg(l,i) ,mapxy(ipg(2,i) ,ipg(3,i))) = rhog(i) 
scrl(img(l,i) ,mapxy(img(2,i) ,img(3,i))) = CONJG [rhog(i)] 

END 

CALL ParallelFFT3D("INV" ,scrl,scr2) 
n(l:NP,l:N y ,l:N z ) = REAL [scr2(l : N£, 1 : N y , 1 : JV,)] 

MODULE FWFFT 

scr2(l:JVP,l:iV y ,l:iV B ) = n(l : N£, 1 : N y , 1 :N Z ) 
CALL ParallelFFT3D("FW" ,scrl,scr2) 
FOR i=l:iVg 

rhog(i) — s cr 1 ( ipg ( 1 , i ) , mapxy ( ipg ( 2 , i ) , ipg ( 3 , i ) ) ) 

END 

Due to the mapping of the y and z direction in Fourier space onto a single dimension, 
input and output array of the parallel Fourier transform do have different shapes. 

MODULE Density 
rto(l:JV x M:JV y ,l:JV z ) = 
FOR i=l:iVb,2 

scrl(l:]V X) l:< c L) = 
FOR j = l:iV£ w 

scrl(ipg(l,i) , mapxy (ipg(2,i) ,ipg(3,i))) = & 

& c(j,i) + I * c(j,i+l) 
scrl(img(l,i) , mapxy (img(2,i) ,img(3,i))) = k 
k C0NJG[c(j,i) + I * c(j,i+l)] 

END 

CALL ParallelFFT3D("INV" ,scrl,scr2) 
rho(l:ArP,l:iVy,l:iV z ) = rho(l : N£, 1 : N y , 1 : N z ) + k 

& REAL [scr2(l: ATP, l:Ar y ,l:AT z )]**2 + k 

k IMAG[scr2(l:7VP,l:AT y ,l:Ar z )]**2 

END 

MODULE VPSI 
FOR i=l:JV b ,2 

scrl(l:JV x ,l:iV p P L) = 
FOR j=l -.N^ 

scrl(ipg(l,i) , mapxy (ipg(2,i) ,ipg(3,i))) = k 

k c(j,i) + I * c(j,i+l) 
scrl(img(l,i) , mapxy (img(2,i) ,img(3,i))) = k 
k C0NJG[c(j,i) + I * c(j,i+l)] 



88 



END 

CALL ParallelFFT3D ( "INV" , scr 1 , scr2) 
scr2(l:NP,l:N y ,l:N z ) = scr2 (1 : N£ , 1 : N y , 1 : N z ) * k 

& vpot(l:JVP,l:JV y> l:AT z ) 
CALL ParallelFFT3D ( "FW" , scr 1 , scr2) 
FDR j = l:^ w 

FP = scrl(ipg(l,i) ,mapxy(ipg(2,i) ,ipg(3,i))) & 

& + scrl(img(l,i) ,mapxy (img(2, i) ,img(3,i))) 

FM = scrl(ipg(l,i) ,mapxy(ipg(2,i) ,ipg(3,i))) & 

& - scrl(img(l,i) ,mapxy(img(2,i) ,img(3,i))) 

fc(j,i) = f(i) * CMPLX [REAL [FP] , IMAG [FM] ] 

fc(j,i+l) = f(i+l) * CMPLX [IMAG [FP] , -REAL [FM] ] 

END 

END 

The parallel Fourier transform routine can be built from a multiple one-dimensional 
Fourier transform and a parallel matrix transpose. As mentioned above, only one 
dimension of the real space grid is distributed in the CPMD code. This allows to 
combine the transforms in y and z direction to a series of two-dimensional trans- 
forms. The handling of the plane waves in Fourier space breaks the symmetry and 
two different transpose routines are needed, depending on the direction. All the 
communication is done in the routine ParallelTranspose. This routine consists 
of a part where the coefficients are gathered into matrix form, the parallel matrix 
transpose, and a final part where the coefficients are put back according to the 
mapping used. 

MODULE ParallelFFT3D(tag,a,b) 
IF (tag == "INV") THEN 
CALL MLTFFTlD(a) 

CALL ParallelTranspose ( " INV" , b , a) 
CALL MLTFFT2D(b) 
ELSE 

CALL MLTFFT2D(b) 

CALL ParallelTranspose ( "FW" , b , a) 
CALL MLTFFTlD(a) 
END IF 

All other parts of the program use the same patterns for the parallelization as the 
ones shown in this section. 

3.9.4 Limitations 

Two types of limitations can be encountered when trying to run a parallel code on a 
computer. Increasing the number of processors working on a problem will no longer 
lead to a faster calculation or the memory available is not sufficient to perform a 
calculation, independently on the number of processors available. The first type of 
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Figure 13. Maximal theoretical speedup for a calculation with a real space grid of dimension 100 
(solid line) . Effective speedup for a 32 water molecule system with an energy cutoff of 70 Rydbcrg 
and a real space grid of dimension 100 (dotted line with diamonds) 



limitation is related to bad load-balancing or the computation becomes dominated 
by the non-scaling part of the communication routines. Load-balancing problems 
in the CPMD code are almost exclusively due to the distribution of the real space 
arrays. Only the x coordinate is distributed. There are typically of the order of 
100 grid points in each direction. Figure 13 shows the maximal theoretical speedup 
for a calculation with a real space grid of dimension 100. The steps are due to the 
load-balancing problems initiated by the granularity of the problem (the dimension 
is an integer value). No further speedup can be achieved once 100 processors are 
reached. The second curve in Fig. 13 shows actual calculations of the full CPMD code. 
It is clearly shown that the load balancing problem in the Fourier transforms affects 
the performance of this special example. Where this steps appear and how severe 
the performance losses are depends of course on the system under consideration. 

To overcome this limitation a method based on processor groups has been im- 
plemented into the code. For the two most important routines where the real space 
grid load balancing problem appears, the calculation of the charge density and the 
application of the local potential, a second level of parallelism is introduced. The 
processors are arranged into a two-dimensional grid and groups are build according 
to the row and column indices. Each processor is a member of its column group 
(colgrp) and its row group (rowgrp). In a first step a data exchange in the column 
group assures that all the data needed to perform Fourier transforms within the 
row groups are available. Then each row group performs the Fourier transforms 
independently and in the end another data exchange in the column groups rebuilds 
the original data distribution. This scheme (shown in the pseudo code for the den- 
sity calculation) needs roughly double the amount of communication. Advantages 
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are the improved load-balancing for the Fourier transforms and the bigger data 
packages in the matrix transposes. The number of plane waves in the row groups 
(JVp^) is calculated as the sum over all local plane waves in the corresponding 
column groups. 

MODULE Density 
rho(l:JVf,l:iV y ,l:]V z ) = 
FOR i=l:iVb,2*Pc 

CALL ParallelTranspose (c ( : ,i) ,colgrp) 

scrl(l:Ar x)1:i vPW ilir ) = 

FOR j = l:A^w 

scrl(ipg(l,i) ,mapxy(ipg(2,i) ,ipg(3,i))) = k 

k c(j,i) + I * c(j,i+l) 
scrl(img(l,i) ,mapxy(img(2,i) ,img(3,i))) = k 

k C0NJG[c(j,i) + I * c(j,i+l)] 

END 

CALL ParallelFFT3D ( " INV" , scr 1 , scr2 , rowgrp) 
rho(l:^P r ,l:^ y ,l:iV z ) = rho (1 : iVf , 1 : N y , 1 : N z ) + k 

k REAL[scr2(l:^|,l:^Vy,l:^Vz)]**2 + k 

k IMAG[scr2(l:A^|,l:^Vy,l:^V z )]**2 

END 

CALL GlobalSum(rho,colgrp) 

The use of two task groups in the example shown in Fig. 13 leads to an increase of 
speedup for 256 processors from 120 to 184 on a Cray T3E/600 computer. 

The effect of the non-scalability of the global communication used in CPMD is 
shown in Fig. 14. This example shows the percentage of time used in the global 
communication routines (global sums and broadcasts) and the time spend in the 
parallel Fourier transforms for a system of 04 silicon atoms with a energy cutoff of 
12 Rydberg. It can clearly be seen that the global sums and broadcasts do not scale 
and therefore become more important the more processors are used. The Fourier 
transforms on the other hand scale nicely for this range of processors. Where 
the communication becomes dominant depends on the size of the system and the 
performance ratio of communication to cpu. 

Finally, the memory available on each processor may become a bottleneck for large 
computations. The replicated data approach for some arrays adapted in the im- 
plementation of the code poses limits on the system size that can be processed on 
a given type of computer. In the outline given in this chapter there are two types 
of arrays that scale quadratically in system size that a replicated. The overlap 
matrix of the projectors with the wavefunctions (f nl) and the overlap matrices of 
the wavefunctions themselves (smat). The fnl matrix is involved in two types of 
calculations where the parallel loop goes either over the bands or the projectors. 
To avoid communication, two copies of the array are kept on each processor. Each 
copy holds the data needed in one of the distribution patterns. This scheme needs 
only a small adaptation of the code described above. 

The distribution of the overlap matrices (smat) causes some more problems. In 
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Figure 14. Percentage of total CPU time spend in global communication routines (solid line) and 
in Fourier transform routines (dashed line) for a system of 64 silicon atoms on a Cray T3E/600 
computer. 

addition to the adaptation of the overlap routine, also the matrix multiply routines 
needed for the orthogonalization step have to be done in parallel. Although there 
arc libraries for these tasks available the complexity of the code is considerably 
increased. 

3.9.5 Summary 

Efficient parallel algorithms for the plane wave-pseudopotential density functional 
theory method exist. Implementations of these algorithms are available and were 
used in most of the large scale applications presented at the end of this paper 
(Sect. 5). Depending on the size of the problem, excellent speedups can be achieved 
even on computers with several hundreds of processors. The limitations presented 
in the last paragraph are of importance for high end applications. Together with 
the extensions presented, existing plane wave codes are well suited also for the next 
generation of supercomputers. 

4 Advanced Techniques: Beyond . . . 

4-1 Introduction 

The discussion up to this point revolved essentially around the "basic" ah initio 
molecular dynamics methodologies. This means in particular that classical nuclei 
evolve in the electronic ground state in the microcanonical ensemble. This com- 
bination allows already a multitude of applications, but many circumstances exist 
where the underlying approximations are unsatisfactory. Among these cases are 
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situations where 



• it is necessary to keep temperature and /or pressure constant (such as during 
journeys in phase diagrams or in the investigation of solid-state phase transi- 
tions), 

• there is a sufficient population of excited electronic states (such as in materials 
with a small or vanishing electronic gap) or dynamical motion occurs in a single 
excited states (such as after photoexcitation events), 

• light nuclei are involved in crucial steps of a process (such as in studies of 
proton transfer or muonium impurities). 

In the following subsections techniques are introduced which transcede these limi- 
tations. Thus, the realm of ah initio molecular dynamics is considerably increased 
beyond the basic setup as discussed in general terms in Sect. 2 and concerning 
its implementation in Sect. 3. The presented "advanced techniques" are selected 
because they are available in the current version of the CPMD package 142 , but their 
implementation is not discussed in detail here. 

4-2 Beyond Microcanonics 
4-2.1 Introduction 

In the framework of statistical mechanics all ensembles can be formally obtained 
from the microcanonical or NVE ensemble - where particle number, volume and 
energy are the external thermodynamic control variables - by suitable Laplace 
transforms of its partition function; note that V is used for volume when it comes 
to labeling the various ensembles in Sect. 4 and its subsections. Thermodynam- 
ically this corresponds to Lcgcndre transforms of the associated thermodynamic 
potentials where intensive and extensive conjugate variables are interchanged. In 
thermodynamics, this task is achieved by a "sufficiently weak" coupling of the 
original system to an appropriate infinitely large bath or reservoir via a link that 
establishes thermodynamic equilibrium. The same basic idea is instrumental in 
generating distribution functions of such ensembles by computer simulation 98 » 250 . 
Here, two important special cases are discussed: thermostats and barostats, which 
are used to impose temperature instead of energy and / or pressure instead of 
volume as external control parameters 12 > 445 > 270 < 585 > 217 . 

4-2.2 Imposing Temperature: Thermostats 

In the limit of ergodic sampling the ensemble created by standard molecular dynam- 
ics is the microcanonical or NVE ensemble where in addition the total momentum 
is conserved 12 > 270 < 217 . Thus, the temperature is not a control variable in the New- 
tonian approach to molecular dynamics and whence it cannot be preselected and 
fixed. But it is evident that also within molecular dynamics the possibility to con- 
trol the average temperature (as obtained from the average kinetic energy of the 
nuclei and the energy equipartition theorem) is welcome for physical reasons. A 
deterministic algorithm of achieving temperature control in the spirit of extended 
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system dynamics 14 by a sort of dynamical friction mechanism was devised by Nose 

and Hoover 442,443,444,307^ gee e g Refe 12,445,270,585,217 for reviews of this well _ 

established technique. Thereby, the canonical or NVT ensemble is generated in 
the case of ergodic dynamics. 

As discussed in depth in Sect. 2.4, the Car-Parrinello approach to ab initio 
molecular dynamics works due to a dynamical separation between the physical 
and fictitious temperatures of the nuclear and electronic subsystems, respectively. 
This separability and thus the associated metastability condition breaks down if the 
electronic excitation gap becomes comparable to the thermal energy or smaller, that 
is in particular for metallic systems. In order to satisfy nevertheless adiabaticity in 
the sense of Car and Parrincllo it was proposed to couple separate thermostats 583 to 
the classical fields that stem from the electronic degrees of freedom 74 > 2 °4 Finally, 
the (long-term) stability of the molecular dynamics propagation can be increased 
due to the same mechanism, which enables one to increase the time step that still 
allows for adiabatic time evolution 638 . Note that these technical reasons to include 
additional thermostats are by construction absent from any Born-Oppcnhcimcr 
molecular dynamics scheme. 

It is well-known that the standard Nose Hoover thermostat method suffers from 
non-ergodicity problems for certain classes of Hamiltonians, such as the harmonic 
oscillator 307 . A closely related technique, the so-called Nose-Hoover-chain ther- 
mostat 388 , cures that problem and assures ergodic sampling of phase space even 
for the pathological harmonic oscillator. This is achieved by thermostatting the 
original thermostat by another thermostat, which in turn is thermostatted and so 
on. In addition to restoring ergodicity even with only a few thermostats in the 
chain, this technique is found to be much more efficient in imposing the desired 
temperature. 

Nose-Hoover-chain thermostatted Car-Parrinello molecular dynamics was in- 
troduced in Ref. 638 . The underlying equations of motion read 

M/ R/ = -V I E KS - M/^R/ (268) 



J2Mi-Rj~gk B T 



Qlik = [Qfe-i4 2 -i - k B T] - Qltkik+i (1 - <W) where k = 2,...,K 
for the nuclear part and 

= -Hf 5 ^ + ^ij4>j - ffli^ (269) 



Qlm = 2 



Qtm = 



J><&|&)-T e ° 

. i 

i 



- Qtmm 



QUvf-i - 



0c 



- QiVlVl+i (1 - Sil) where I = 2, . . . , L 



for the electronic contribution. These equations are written down in density func- 
tional language (see Eq. (75) and Eq. (81) for the definitions of E KS and H^ s , 
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respectively), but completely analogues expressions are operational if other elec- 
tronic structure approaches are used instead. Using separate thermostatting baths 
and {r/i}, chains composed of K and L coupled thermostats are attached to 
the nuclear and electronic equations of motion, respectively. 

By inspection of Eq. (268) it becomes intuitively clear how the thermostat works: 
£1 can be considered as a dynamical friction coefficient. The resulting "dissipative 
dynamics" leads to non-Hamiltonian flow, but the friction term can aquire positive 
or negative sign according to its equation of motion. This leads to damping or 
acceleration of the nuclei and thus to cooling or heating if the instantaneous kinetic 
energy of the nuclei is higher or lower than k^T which is preset. As a result, 
this extended system dynamics can be shown to produce a canonical ensemble 
in the subspace of the nuclear coordinates and momenta. In spite of being non- 
Hamiltonian, Nose-Hoover (-chain) dynamics is also distinguished by conserving 
an energy quantity of the extended system, see Eq. (272). 

The desired average physical temperature is given by T and g denotes the num- 
ber of dynamical degrees of freedom to which the nuclear thermostat chain is cou- 
pled (i.e. constraints imposed on the nuclei have to be subtracted). Similarly, T ° is 
the desired fictitious kinetic energy of the electrons and 1 / (3 C is the associated tem- 
perature. In principle, /3 e should be chosen such that l//3 e = 2T®/N e where N e is 
the number of dynamical degrees of freedom needed to parameterize the wavefunc- 
tion minus the number of constraint conditions. It is found that this choice requires 
a very accurate integration of the resulting equations of motion (for instance by us- 
ing a high-order Suzuki- Yoshida integrator, see Sect. VI. A in Ref. 638 ). However, 
relevant quantities are rather insensitive to the particular value so that N c can be 
replaced heuristically by which is the number of orbitals fa used to expand the 
wavefunction 638 . 

The choice of the "mass parameters" assigned to the thermostat degrees of 
freedom should be made such that the overlap of their power spectra and the ones 
the thermostatted subsystems is maximal 74 ' 638 . The relations 

n n 

« = Q °=ik (271) 

assures this if uj u is a typical phonon or vibrational frequency of the nuclear subsys- 
tem (say of the order of 2000 to 4000 cm -1 ) and uj c is sufficiently large compared 
to the maximum frequency w™ ax of the nuclear power spectrum (say 10 000 cm -1 
or larger). The integration of these equations of motion is discussed in detail in 
Ref. 638 using the velocity Verlet / rattle algorithm. 

In some instances, for example during equilibration runs, it is advantageous to 
go one step further and to actually couple one chain of Nose-Hoover thermostats 
to every individual nuclear degree of freedom akin to what is done in path integral 
molecular dynamics simulations 637,644,646^ gee a ^ g0 g ec ^ 4 4 xhis so-called "mas- 
sive thermostatting approach" is found to accelerate considerably the expensive 
equilibration periods within ab initio molecular dynamics, which is useful for both 
Car-Parrinello and Born-Oppcnhcimcr dynamics. 
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In classical molecular dynamics two quantities are conserved during a simula- 
tion, the total energy and the total momentum. The same constants of motion 
apply to (exact) microcanonical Born Oppenheimcr molecular dynamics because 
the only dynamical variables are the nuclear positions and momenta as in classi- 
cal molecular dynamics. In microcanonical Car-Parrincllo molecular dynamics the 
total energy of the extended dynamical system composed of nuclear and electronic 
positions and momenta, that is E cons as defined in Eq. (48), is also conserved, see 
e.g. Fig. 3 in Sect. 2.4. There is also a conserved energy quantity in the case of ther- 
mostatted molecular dynamics according to Eq. (268)-(269). Instead of Eq. (48) 
this constant of motion reads 



eZI = E" |* ) + E + £ KS [{*}, {*/}] 

i I 

+E^?+Ej + 2T eV 

1=1 1 1=2 Pe 

K K 

+ Eo^ + E fc Bra+5fcBTa (272) 



2 

k=l k=2 

for Nose-Hoover-chain thermostatted canonical Car-Parrinello molecular dynam- 
ics 638 . 

In microcanonical Car-Parrinello molecular dynamics the total nuclear momen- 
tum P n is no more a constant of motion as a result of the fictitious dynamics of the 
wavefunction; this quantity as well as other symmetries and associated invariants 
are discussed in Ref. 467 . However, a generalized linear momentum which embraces 
the electronic degrees of freedom 

occ 

P C P = Pn + Pc = E P ' + E » - Vr 0i) + C.C. (273) 



can be defined 467 > 436 ; p 7 = MjRj. This quantity is a constant of motion in 
unthermostatted Car-Parrinello molecular dynamics due to an exact cancellation of 
the nuclear and electronic contributions 467 > 436 . As a result, the nuclear momentum 
P n fluctuates during such a run, but in practice P n is conserved on the average as 
shown in Fig. 1 of Ref. 436 . This is analogues to the behavior of the physical total 
energy E p h ys Eq. (49) , which fluctuates slightly due to the presence of the fictitious 
kinetic energy of the electrons T e Eq. (51). 

As recently outlined in detail it is clear that the coupling of more than one 
thermostat to a dynamical system, such as done in Eq. (268)-(269), destroys the 
conservation of momentum 436 , i.e. Pep is no more an invariant. In unfavorable 
cases, in particular in small-gap or metallic regimes where there is a substantial 
coupling of the nuclear and electronic subsystems, momentum can be transferred 
to the nuclear subsystem such that P n grows in the course of a simulation. This 
problem can be cured by controlling the nuclear momentum (using e.g. scaling or 
constraint methods) so that the total nuclear momentum P n remains small 436 . 
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4-2.3 Imposing Pressure: Barostats 

Keeping the pressure constant is a desirable feature for many applications of molec- 
ular dynamics. The concept of barostats and thus constant-pressure molecular dy- 
namics was introduced in the framework of extended system dynamics by Hans An- 
dersen 14 , see e.g. Refs. 12 > 270 < 585 > 217 f or introductions. This method was devised to 
allow for isotropic fluctuations in the volume of the supercell. A powerful extension 
consists in also allowing for changes of the shape of the supercell to occur as a result 
of applying external pressure 459 > 46 °. 4 6i,678^ i ncnlc r m g the possibility of non-isotropic 
external stress 460 ; the additional fictitious degrees of freedom in the Parrinello- 
Rahman approach 459 > 460 > 4 6i are the lattice vectors of the supercell, whereas the 
strain tensor is the dynamical variable in the Wentzcovitch approach 678 . These 
variable-cell approaches make it possible to study dynamically structural phase 
transitions in solids at finite temperatures. With the birth of ah initio molecu- 
lar dynamics both approaches were combined starting out with isotropic volume 
fluctuations 94 a la Andersen 14 and followed by Born-Oppenheimer 681 > 682 and 
Car-Parrinello 201 > 202 » 55 » 56 variable-cell techniques. 

The basic idea to allow for changes in the cell shape consists in constructing 
an extended Lagrangian where the primitive Bravais lattice vectors ai, a 2 and a 3 
of the simulation cell are additional dynamical variables similar to the thermostat 
degree of freedom £, see Eq. (268). Using the 3x3 matrix h = [ai,a 2 ,a 3 ] (which 
fully defines the cell with volume Q) the real-space position R/ of a particle in this 
original cell can be expressed as 

R/ = hSj (274) 

where S/ is a scaled coordinate with components S/.„ e [0, 1] that defines the 
position of the 7th particle in a unit cube (i.e. fl un it = 1) which is the scaled 
cell 459 ' 460 ) see Sect. 3.1 for some definitions. The resulting metric tensor Q = h*h 
converts distances measured in scaled coordinates to distances as given by the 
original coordinates according to Eq. (106) and periodic boundary conditions are 
applied using Eq. (107). 

In the case of ah initio molecular dynamics the orbitals have to be expressed 
suitably in the scaled coordinates s = h~ 1 r. The normalized original orbitals <f>i(r) 
as defined in the unsealed cell h are transformed according to 

&(r) = ^=&(s) (275) 

satisfying 

f drti(r)Mr)= [ <M?(s)&(s) (276) 
so that the resulting charge density is given by 

n(r) = ±n(s) . (277) 

in the scaled cell, i.e. the unit cube. Importantly, the scaled fields <fii(s) and thus 
their charge density n(s) do not depend on the dynamical variables associated to 
the cell degrees of freedom and thus can be varied independently from the cell; the 
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original unsealed fields 4>i{v) do depend on the cell variables h via the normalization 
by the cell volume = det h as evidenced by Eq. (275). 

After these preliminaries a variable-cell extended Lagrangian for ab initio molec- 
ular dynamics can be postulated 202 > 201 > 55 

C = J2 » (<M S ) |<M S ) ) - E KS [{&}, {hS/}] 

i 

ij 

+ S \ Ml (S/SSj) + \w Tr h'h - p O , (278) 

with additional nine dynamical degrees of freedom that are associated to the lat- 
tice vectors of the supercell h. This constant-pressure Lagrangian reduces to the 
constant-volume Car-Parrinello Lagrangian, see e.g. Eq. (41) or Eq. (58), in the 
limit h — > of a rigid cell (apart from a constant term p ft). Here, p defines the 
externally applied hydrostatic pressure, W defines the fictitious mass or inertia pa- 
rameter that controls the time-scale of the motion of the cell h and the interaction 
energy E KS is of the form that is defined in Eq. (75). In particular, this Lagrangian 
allows for symmetry-breaking fluctuations - which might be necessary to drive 
a solid-state phase transformation - to take place spontaneously The resulting 
equations of motion read 

Ml S IiU = -£__ (h*)" 1 - Mr £ ]T G-Hfc, A. (279) 

v—1 ' v— 1 8=1 

5E KS 

M<Ms) = +I>«&(s) ( 28 °) 

3 

WKv = nY,{K?-pS™)(P)^ , (281) 

8=1 

where the total internal stress tensor 

K1 = E Ml (S/^S/) us + n - ( 282 ) 

is the sum of the thermal contribution due to nuclear motion at finite temperature 
and the electronic stress tensor 440 > 441 n which is defined in Eq. (189) and the 
following equations, see Sect. 3.4. 

Similar to the thermostat case discussed in the previous section one can rec- 
ognize a sort of frictional feedback mechanism. The average internal pressure 
((1/3) Tr II tot ) equals the externally applied pressure p as a result of maintain- 
ing dynamically a balance between p S and the instantaneous internal stress II tot 
by virtue of the friction coefficient oc Q in Eq. (279). Ergodic trajectories obtained 
from solving the associated ab initio equations of motion Eq. (279)-(281) lead to 
a sampling according to the isobaric-isoenthalpic or NpH ensemble. However, the 
generated dynamics is fictitious similar to the constant-temperature case discussed 
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in the previous section. The isobaric-isothermal or NpT ensemble is obtained by 
combining barostats and thermostats, see Ref. 389 for a general formulation and 
Rcf. 391 for reversible integration schemes. 

An important practical issue in isobaric ab initio molecular dynamics simula- 
tions is related to problems caused by using a finite basis set, i.e. "incomplctc- 
basis-set" or Pulay-type contributions to the stress, see also Sect. 2.5. Using a 
finite plane wave basis (together with a finite number of k-points) in the presence 
of a fluctuating cell 245 > 211 one can either fix the number of plane waves or fix the 
energy cutoff; see Eq. (122) for their relation according to a rule-of-thumb. A 
constant number of plane waves implies no Pulay stress but a decreasing precision 
of the calculation as the volume of the supercell increases, whence leading to a sys- 
tematically biased (but smooth) equation of state. The constant cutoff procedure 
has better convergence properties towards the infinite-basis-set limit 245 . However, 
it produces in general unphysical discontinuities in the total energy and thus in the 
equation of state at volumes where the number of plane waves changes abruptly, 
see e.g. Fig. 5 in Ref. 211 . 

Computationally, the number of plane waves has to be fixed in the framework 
of Car-Parrinello variable-cell molecular dynamics 94 > 202 > 201 > 55 ; whereas the energy 
cutoff can easily be kept constant in Born-Oppcnhcimer approaches to variable- 
cell molecular dynamics 681 > 682 . Sticking to the Car-Parrinello technique a practical 
remedy 202 > 55 to this problem consists in modifying the electronic kinetic energy 
term Eq. (173) in a plane wave expansion Eq. (172) of the Kohn-Sham functional 
E KS Eq. (75) 

£kin = l>E^ G | 2 |Q(q)| 2 , (283) 

i q 

where the unsealed G and scaled q = 27rg reciprocal lattice vectors arc interrelated 
via the cell h according to Eq. (Ill) (thus Gr = qs) and the cutoff Eq. (121) 
is defined as (1/2) |G| 2 < E cut for a fixed number of q- vectors, see Sect. 3.1. 
The modified kinetic energy at the T-point of the Brillouin zone associated to the 
supercell reads 

^ki„ = E/*EJ|g (A^c°ut) 2 Ic(q)| 2 (284) 



2 
q 



= |G| 2 + ^ 1 + erf 



^|G| 2 -£ c ° u ff t 



(285) 



where A, a and E^ t are positive definite constants and the number of scaled vectors 
q, that is the number of plane waves, is strictly kept fixed. 

In the limit of a vanishing smoothing (A — > 0; a — > oo) the constant number of 
plane wave result is recovered. In limit of a sharp step function (A — > oo; a — > 0) 
all plane waves with (1/2) |G| 2 >• E^ t have a negligible weight in E^ n and are 
thus effectively suppressed. This situation mimics a constant cutoff calculation at 
an "effective cutoff" of w -E£ut within a constant number of plane wave scheme. For 
this trick to work note that E cut 3> £^ut nas to be satisfied. In the case ^4)0 the 
electronic stress tensor II given by Eq. (189) features an additional term (due to 
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changes in the "effective basis set" as a result of variations of the supercell), which 
is related to the Pulay stress 219 < 660 . 

Finally, the strength of the smoothing A)0 should be kept as modest as possible 
since the modification Eq. (284) of the kinetic energy leads to an increase of the 
highest frequency in the electronic power spectrum oc A. This implies a decrease 
of the permissible molecular dynamics time step Ai max according to Eq. (55). It 
is found that a suitably tuned set of the additional parameters (A, a, E^ t ) leads 
to an efficiently converging constant-pressure scheme in conjunction with a fairly 
small number of plane waves 202 > 55 . Note that the cutoff was kept strictly con- 
stant in applications of the Born-Oppenheimer implementation 679 of variable-cell 
molecular dynamics 681 > 682 ; but the smoothing scheme presented here could be 
implemented in this case as well. An efficient method to correct for the discontinu- 
ities of static total energy calculations performed at constant cutoff was proposed 
in Ref. 211 . Evidently, the best way to deal with the incomplete-basis-set problem 
is to increase the cutoff such that the resulting artifacts become negligible on the 
physically relevant energy scale. 

4-3 Beyond Ground States 
4- 3.1 Introduction 

Extending ab initio molecular dynamics to a single non-interacting excited state is 
straightforward in the framework of wavefunction-based methods such as Hartree- 

Fock 365,254,191,379,281,284,316,293 ; generalized valence bond (GVB) 2 82,283,228,229,230 ; 

complete active space SCF (CASSCF) 566 < 567 ; or full configuration interaction 
(FCI) 372 approaches, see Sect. 2.7. However, these methods are computationally 
quite demanding - given present-day algorithms and hardware. Promising steps 
in the direction of including several excited states and non-adiabatic couplings are 

also made 385,386,387,71^ 

Density functional theory offers an alternative route to approximately solving 
electronic structure problems and recent approaches to excited-state properties 
within this framework look promising. In the following, two limiting and thus 
certainly idealistic situations will be considered, which are characterized by either 

• many closely-spaced excited electronic states with a broad thermal distribution 
of fractional occupation numbers, or by 

• a single electronic state that is completely decoupled from all other states. 

The first situation is encountered for metallic systems with collective excitations or 
for materials at high temperatures compared to the Fermi temperature. It is noted 
in passing that associating fractional occupation numbers to one-particle orbitals is 
also one route to go beyond a single-determinant ansatz for constructing the charge 
density 458 > 168 . The second case applies for instance to large-gap molecular systems 
which complete a chemical reaction in a single excited state as a result of e.g. a 
vertical homo / lumo or instantaneous one-particle / one-hole photoexcitation. 
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4-3.2 Many Excited States: Free Energy Functional 

The free energy functional approach to excited-state molecular dynamics 5,7 is a 
mean field approach similar in spirit to Ehrenfest molecular dynamics, see Sect. 2.2. 
The total wavcfunction is first factorized into a nuclear and an electronic wave- 
function Eq. (3) followed by taking the classical limit for the nuclear subsystem. 
Thus, classical nuclei move in the average field as obtained from averaging over all 
electronic states Eq. (25). A difference is that according to Ehrenfest molecular 
dynamics the electrons are propagated in real time and can perform non-adiabatic 
transitions by virtue of direct coupling terms oc d kl between all states \tfc subject 
to energy conservation, see Sect. 2.2 and in particular Eqs. (27)-(29). The average 
force or Ehrenfest force is obtained by weighting the different states k according to 
their diagonal density matrix elements (that is oc |cfc(i)| 2 in Eq. (27)) whereas the 
coherent transitions are driven by the off-diagonal contributions (which are oc c* k ci, 
see Eq. (27)). 

In the free energy approach 5 ' 7 , the excited states are populated according to 
the Fermi Dirac (finite temperature equilibrium) distribution which is based on 
the assumption that the electrons "equilibrate" more rapidly than the timescale 
of the nuclear motion. This means that the set of electronic states evolves at a 
given temperature "isothermally" (rather than adiabatically) under the inclusion 
of incoherent electronic transitions at the nuclei move. Thus, instead of comput- 
ing the force acting on the nuclei from the electronic ground-state energy it is 
obtained from the electronic free energy as defined in the canonical ensemble. By 
allowing such electronic transitions to occur the free energy approach transcends 
the usual Born-Oppenheimer approximation. However, the approximation of an 
instantaneous equilibration of the electronic subsystem implies that the electronic 
structure at a given nuclear configuration {R/} is completely independent from pre- 
vious configurations along a molecular dynamics trajectory. Due to this assumption 
the notion "free energy Born-Oppenheimer approximation" was coined in Ref. 101 
in a similar context. Certain non-equilibrium situations can also be modeled within 
the free energy approach by starting off with an initial orbital occupation pattern 
that does not correspond to any temperature in its thermodynamic meaning, see 
e.g. Refs. 57 °. 572 > 574 for such applications. 

The free energy functional as defined in Ref. 5 is introduced most elegantly 7 ' 9 
by starting the discussion for the special case of non-interacting Fermions 

". = 4 V2 -S^ (286) 

in a fixed external potential due to a collection of nuclei at positions {R/}. The 
associated grand partition function and its thermodynamic potential ("grand free 
energy") are given by 

3 B (/iVT) = det 2 (1 + exp [-0 (H s - fj,)]) (287) 

fl s {tJ,VT) = -k B TlnZ s {nVT) , (288) 
where \x is the chemical potential acting on the electrons and the square of the 
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determinant stems from considering the spin unpolarized special case only. This 
reduces to the well-known grand potential expression 

n s (fiVT) = -2fc B Tlndet(l + exp [-0{H a - //)]) 

= -2fc B T^ln(l + exp -/? (e^ - fi) ) (289) 

i 

for non-interacting spin-1/2 Fermions where {e^} are the eigenvalues of a one 
particle Hamiltonian such as Eq. (286); here the standard identity IndetM = 
TrlnM was invoked for positive definite M. 

According to thermodynamics the Helmholtz free energy T(NVT) associated 
to Eq. (288) can be obtained from an appropriate Legendre transformation of the 
grand free energy Q(fiVT) 

T S {NVT) = n s (fiVT) + fiN + J2 ]vt Z ' Z '' | (290) 

by fixing the average number of electrons N and determining \x from the conven- 
tional thermodynamic condition 



VT 



In addition, the internuclcar Coulomb interactions between the classical nuclei were 
included at this stage in Eq. (290). Thus, derivatives of the free energy Eq. (290) 
with respect to ionic positions — V/JF S define forces on the nuclei that could be used 
in a (hypothetical) molecular dynamics scheme using non-interacting electrons. 

The interactions between the electrons can be "switched on" by resorting to 
Kohn-Sham density functional theory and the concept of a non-interacting refer- 
ence system. Thus, instead of using the simple one-particle Hamiltonian Eq. (286) 
the effective Kohn-Sham Hamiltonian Eq. (83) has to be utilized. As a result, the 
grand free energy Eq. (287) can be written as 

n KS (nVT) = -2fc B Tln [det (l + exp [-/? (H KS - p)] )] (292) 

" KS H 72 -£jidV v " (r, + W (293) 

H KS <f>i = (294) 

where f2 xc is the exchange-correlation functional at finite temperature. By virtue of 
Eq. (289) one can immediately sec that Q KS is nothing else than the "Fermi-Dirac 
weighted sum" of the bare Kohn-Sham eigenvalues {ti}. Whence, this term is the 
extension to finite temperatures of the "band-structure energy" (or of the "sum 
of orbital energies" in the analogues Hartree-Fock case 604 > 418 ) contribution to the 
total electronic energy, see Eq. (86). 

In order to obtain the correct total electronic free energy of the interacting 
electrons as defined in Eq. (86) the corresponding extra terms (properly generalized 
to finite temperatures) have to be included in Q KS . This finally allows one to write 
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down the generalization of the Hclmholtz free energy of the interacting many- 
electron case 



\ J dr Vu( 



i (J 

.■!«lvH o xi - I dr 5 -^yn(r) (295) 

in the framework of a Kohn-Sham-like formulation. The corresponding one- 
particle density at the T-point is given by 

(296) 



i(r) = £ MP) \Mr) 



/<(/?) = (1 + exp [/?(*-/*)]) 1 , (297) 

where the fractional occupation numbers {/*} are obtained from the Fermi-Dirac 
distribution at temperature T in terms of the Kohn Sham eigenvalues {ei}. Finally, 
ab initio forces can be obtained as usual from the nuclear gradient of JF KS , which 
makes molecular dynamics possible. 

By construction, the total free energy Eq. (295) reduces to that of the non- 
interacting toy model Eq. (290) once the electron electron interaction is switched 
off. Another useful limit is the ground-state limit [3 — > oo where the free energy 
F KS (NVT) yields the standard Kohn-Sham total energy expression E KS as de- 
fined in Eq. (86) after invoking the appropriate limit f2 xc — * E xc as T — > 0. Most 
importantly, stability analysis 5 ' 7 of Eq. (295) shows that this functional shares the 
same stationary point as the exact finite-temperature functional due to Mermin 424 , 
see e.g. the textbooks 458 ' 168 for introductions to density functional formalisms at 
finite temperatures. This implies that the self-consistent density, which defines 
the stationary point of !F KS , is identical to the exact one. This analysis reveals 
furthermore that, unfortunately, this stationary point is not an extremum but a 
saddle point so that no variational principle and, numerically speaking, no direct 
minimization algorithms can be applied. For the same reason a Car-Parrinello 
fictitious dynamics approach to molecular dynamics is not a straightforward op- 
tion, whereas Born-Oppenheimcr dynamics based on diagonalization can be used 
directly. 

The band-structure energy term is evaluated in the CPMD package 142 by di- 
agonalizing the Kohn-Sham Hamiltonian after a suitable "preconditioning" 5 , see 
Sect. 3.6.2. Specifically, a second-order Trotter approximation is used 

Tr exp [-0H KS ] = ]T exp [-fa] = ]T Pu (/?) (298) 
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in order to compute first the diagonal elements pa (At) of the "high-temperature" 
Boltzmann operator p(Ar); here At = 0/P and P is the Trotter "time slice" as 
introduced in Sect. 4.4.2. To this end, the kinetic and potential energies can be 
conveniently evaluated in reciprocal and real space, respectively, by using the split 
operator / FFT technique 183 . The Kohn-Sham eigenvalues ej are finally obtained 
from the density matrix via ej = — (1/Ar) \npn(AT). They are used in order to 
compute the occupation numbers {/,}, the density n(r), the band-structure energy 
il KS , and thus the free energy Eq. (295). 

In practice a diagonalization / density-mixing scheme is employed in order to 
compute the self-consistent density n(r). Grossly speaking a suitably constructed 
trial input density n- ln (see e.g. the Appendix of Ref. 571 for such a method) is used 
in order to compute the potential y KS [n; n ]. Then the lowest-order approximant 
to the Boltzmann operator Eq. (300) is diagonalized using an iterative Lanczos- 
type method. This yields an output density n out and the corresponding free energy 

^KS 

[n out ]. Finally, the densities are mixed and the former steps are iterated until a 
stationary solution n sc f of .F KS [n sc f] is achieved, see Sect. 3.6.4 for some details on 
such methods. Of course the most time-consuming part of the calculation is in the 
iterative diagonalization. In principle this is not required, and it should be possible 
to compute the output density directly from the Fermi-Dirac density matrix even 
in a linear scaling scheme 243 , thus circumventing the explicit calculation of the 
Kohn-Sham eigenstates. However, to date efforts in this direction have failed, or 
given methods which are too slow to be useful 9 . 

As a method, molecular dynamics with the free energy functional is most ap- 
propriate to use when the excitation gap is either small, or in cases where the 
gap might close during a chemical transformation. In the latter case no instabil- 
ities are encountered with this approach, which is not true for ground-state ah 
initio molecular dynamics methods. The price to pay is the quite demanding it- 
erative computation of well-converged forces. Besides allowing such applications 
with physically relevant excitations this method can also be straightforwardly com- 
bined with k-point sampling and applied to metals at "zero" temperature. In 
this case, the electronic "temperature" is only used as a smearing parameter of 
the Fermi edge by introducing fractional occupation numbers, which is known to 
improve greatly the convergence of these ground-state electronic structure calcula- 

t j ong 220,232,185,676,680,343,260,344,414,243 

Finite-temperature expressions for the exchange-correlation functional fi xc are 
available in the literature. However, for most temperatures of interest the correc- 
tions to the ground-state expression are small and it seems justified to use one of 
the various well-established parameterizations of the exchange-correlation energy 
E xc at zero temperature, see Sect. 2.7. 

4-3.3 A Single Excited State: Si-Dynamics 

For large-gap systems with well separated electronic states it might be desirable 
to single out a particular state in order to allow the nuclei to move on the asso- 
ciated excited state potential energy surface. Approaches that rely on fractional 
occupation numbers such as ensemble density functional theories including the 
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free energy functional discussed in the previous section - are difficult to adapt for 
cases where the symmetry and / or spin of the electronic state should be fixed 168 . 
An early approach in order to select a particular excited state was based on intro- 
ducing a "quadratic restoring potential" which vanishes only at the eigenvalue of 
the particular state 417 ' m . 

A method that combines Roothaan's symmetry-adapted wavefunctions with 
Kohn-Sham density functional theory was proposed in Rcf. 214 and used to simulate 
a photoisomerization via molecular dynamics. Viewed from Kohn-Sham theory this 
approach consists in building up the spin density of an open-shell system based on 
a symmetry-adapted wavefunction that is constructed from spin-restricted deter- 
minants (the "microstates"). Viewed from the restricted open-shell Hartree-Fock 
theory d la Roothaan it amounts essentially to replacing Hartree-Fock exchange by 
an approximate exchange-correlation density functional. This procedure leads to 
an orbital-dependent density functional which was formulated explicitely for the 
first-excited singlet state 5i in Ref. 214 . The relation of this approach to previ- 
ous theories is discussed in some detail in Ref. 214 . In particular, the success of the 
closely-related Ziegler-Rauk-Baerends "sum methods" 704 > 150 > 600 wa s an important 
stimulus. More recently several papers 252 < 439 > 193 > 195 . 19 6 appeared that are similar 
in spirit to the method of Ref. 214 . The approach of Refs. 193 > 195 > 196 can be viewed 
as a generalization of the special case (Si state) worked out in Ref. 214 to arbitrary 
spin states. In addition, the generalized method 193 > 195 > 196 was derived within the 
framework of density functional theory, whereas the wavefunction perspective was 
the starting point in Ref. 214 . 

In the following, the method is outlined with the focus to perform molecular 
dynamics in the Si state. Promoting one electron from the homo to the LUMO in a 
closed-shell system with 2n electrons assigned to n doubly occupied orbitals (that is 
spin-restricted orbitals that have the same spatial part for both spin up a and spin 
down j3 electrons) leads to four different excited wavefunctions or determinants, see 
Fig. 15 for a sketch. Two states \t\) and ^2) are energetically degenerate triplets 
t whereas the two states \m\) and \rri2} are not eigenfunctions of the total spin 
operator and thus degenerate mixed states m ("spin contamination"). Note in 
particular that the m states do not correspond - as is well known - to singlet states 
despite the suggestive occupation pattern in Fig. 15. 

However, suitable Clebsch-Gordon projections of the mixed states \m\) and 
|m 2 ) yield another triplet state ^3) and the desired first excited singlet or S\ state 
\s\). Here, the ansatz 214 for the total energy of the S\ state is given by 

E Sl [{&}] = 2E™ [{&}] - E™ [{&}] (301) 

where the energies of the mixed and triplet determinants 

El s [{^}]=T s [n] + J dv Kxt(r)n(r) + \ J dr V H (r)n(r) + E xc [n^, r&] (302) 

E? s {{&}} =T s [n] + J dr V cxt (r)n(r) + \ J dr Vfe(r)n(r) + E xc [n?, nf] (303) 

are expressed in terms of (restricted) Kohn-Sham spin-density functional con- 
structed from the set {&}, cf. Eq. (75). The associated Si wavefunction is given 
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Figure 15. Four possible determinants |ti), 1*2}, |mi) and |rri2} as a result of the promotion of a 
single electron from the HOMO to the LUMO of a closed shell system, see text for further details. 
Taken from Ref. 214 . 



by 

| ax [{&}] ) = y/2\m [{&}] } - 1 1 [{fa}] ) (304) 

where the "microstates" m and £ are both constructed from the same set {4>i} of 
n+1 spin-restricted orbitals. Using this particular set of orbitals the total density 

n(r) = n° (r) + n£ (r) = < (r) + nf (r) (305) 

is of course identical for both the m and t determinants whereas their spin den- 
sities clearly differ, see Fig. 16. Thus, the decisive difference between the m and 
t functionals Eq. (302) and Eq. (303), respectively, comes exclusively from the 
exchange-correlation functional E xc , whereas kinetic, external and Hartree energy 
are identical by construction. Note that this basic philosophy can be generalized 
to other spin-states by adapting suitably the microstates and the corresponding 
coefficients in Eq. (301) and Eq. (304). 

Having defined a density functional for the first excited singlet state the 
corresponding Kohn Sham equations are obtained by varying Eq. (301) using 
Eq. (302) and Eq. (303) subject to the orthonormality constraint Y^ij=i ^-ijd&i I 
4>j) — Sij). Following this procedure the equation for the doubly occupied orbitals 
i = 1 , . . . , n — 1 reads 
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Figure 16. Four patterns of spin densities nf, , n^, and n„ corresponding to the two spin- 
restricted determinants \t) and \m) sketched in Fig. 15, see text for further details. Taken from 
Ref. 214 . 



whereas 



and 




are two different equations for the two singly-occupied open-shell orbitals a and 
b, respectively, see Fig. 15. Note that these Kohn-Sham-like equations fea- 
ture an orbital-dependent exchange-correlation potential where V^"[n^,n^J = 
SExdn^, n^/Sn^ and analogues definitions hold for the (3 and t cases. 

The set of equations Eq. (306)-(308) could be solved by diagonalization of the 
corresponding "restricted open shell Kohn-Sham Hamiltonian" or alternatively by 
direct minimization of the associated total energy functional. The algorithm pro- 
posed in Ref. 240 , which allows to properly and efficiently minimize such orbital- 
dependent functionals including the orthonormality constraints, was implemented 
in the CPMD package 142 . Based on this minimization technique Born-Oppcnheimcr 
molecular dynamics simulations can be performed in the first excited singlet state. 
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The alternative Car-Parrincllo formulation seems inconvenient because the singly 
and doubly occupied orbitals would have to be constrained not to mix. 

4-4 Beyond Classical Nuclei 
4-4-1 Introduction 

Up to this point the nuclei were approximated as classical point particles as custom- 
arily done in standard molecular dynamics. There are, however, many situations 
where quantum dispersion broadening and tunneling effects play an important role 
and cannot be neglected if the simulation aims at being realistic - which is the 
generic goal of ah initio simulations. The ah initio path integral technique 395 and 
its extension to quasiclassical time evolution 411 is able to cope with such situations 
at finite temperatures. It is also implemented in the CPMD package 142 . The central 
idea is to quantize the nuclei using Feynman's path integrals and at the same time 
to include the electronic degrees of freedom akin to ah initio molecular dynamics - 
that is "on-the-fly". The main ingredients and approximations underlying the ah 
initio path integral approach 395,399,644,404 are 

• the adiabatic separation of electrons and nuclei where the electrons are kept in 
their ground state without any coupling to electronically excited states (Born- 
Oppenheimer or "clamped-nuclei" approximation), 

• using a particular approximate electronic structure theory in order to calculate 
the interactions, 

• approximating the continuous path integral for the nuclei by a finite discretiza- 
tion (Trotter factorization) and neglecting the indistinguishability of identical 
nuclei (Boltzmann statistics), and 

• using finite supercells with periodic boundary conditions and finite sampling 
times (finite-size and finite-time effects) as usual. 

Thus, quantum effects such as zero-point motion and tunneling as well as thermal 
fluctuations are included at some preset temperature without further simplifications 
consisting e.g. in quasiclassical or quasiharmonic approximations, restricting the 
Hilbert space, or in artificially reducing the dimensionality of the problem. 

4-4-2 Ab Initio Path Integrals: Statics 

For the purpose of introducing ab initio path integrals 395 it is convenient to start 
directly with Feynman's formulation of quantum-statistical mechanics in terms 
of path integrals as opposed to Schrodinger's formulation in terms of wavefunc- 
tions which was used in Sect. 2.1 in order to derive ab initio molecular dynamics. 
For a general introduction to path integrals the reader is referred to standard text- 
books 187 488,334^ wnereas their use in numerical simulations is discussed for instance 
in Refs 233 . 126 < 542 > 120 > 12 4, 646,407 

The derivation of the expressions for ab initio path integrals is based on assuming 
the non-rclativistic standard Hamiltonian, sec Eq. (2). The corresponding canoni- 
cal partition function of a collection of interacting nuclei with positions R = {R/} 
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and electrons r = {r^} can be written as a path integral 



H = j> VR j> Vrcxp -IJ'ZtCe ({R/(r)}, {Rj(r)}; {^(t)}, {^(r)} 



(309) 



where 



C E = T(R) + V(R) + T(f ) + V(r) + V(R, r) 

"V2 Mj {-*r) + ir J W7^7\ 



Hi,-, 



d7 



y — 



e 2 Zj 



(310) 



denotes the Euclidean Lagrangian. The primes in Eq. (309) indicate that the proper 
sums over all permutations corresponding to Bose-Einstein and/or Fermi-Dirac 
statistics have to be included. It is important to note that in Eq. (309) and (310) 
the positions {R/} and {r^} are not operators but simply functions of the imaginary 
time r e [0, K[3] which parameterizes fluctuations around the classical path. This 
implies that the dots denote here derivatives with respect to imaginary time r as 
defined in Eq. (310). According to Eq. (309) exact quantum mechanics at finite 
temperature T — 1/&B/3 is recovered if all closed paths [{R/}; {r^}] of "length" h(3 
are summed up and weighted with the exponential of the Euclidean action measured 
in units of fi; atomic units will be used again from here on. The partial trace over 
the electronic subsystem can formally be written down exactly 



Z = 



<j> PRcxp - J oIt (T(R) + V(R)) 



with the aid of the influence functional 187 > 334 



■ZfRl 



] = j> Vr cxp - dr (T(r) + V(r) + V(R, r)) 



(311) 



(312) 



Note that Z[H] is a complicated and unknown functional for a given nuclear path 
configuration [{R/}]. As a consequence the interactions between the nuclei become 
highly nonlocal in imaginary time due to memory effects. 

In the standard Born-Oppcnhcimcr or "clamped nuclei" approximation, see 
Rcf. 340 for instance, the nuclei are frozen in some configuration and the complete 
electronic problem is solved for this single static configuration. In addition to the 
nondiagonal correction terms that are already neglected in the adiabatic approxi- 
mation, the diagonal terms are now neglected as well. Thus the potential for the 
nuclear motion is simply defined as the bare electronic eigenvalues obtained from a 
series of fixed nuclear configurations. 

In the statistical mechanics formulation of the problem Eq. (311)— (312) the 
Born-Oppcnhcimcr approximation amounts to a "quenched average" : at imaginary 
time r the nuclei are frozen at a particular configuration R(t) and the electrons 
explore their configuration space subject only to that single configuration. This 
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implies that the electronic degrees of freedom at different imaginary times r and r' 
become completely decoupled. Thus, the influence functional Z [R] has to be local in 
t and becomes particularly simple; a discussion of adiabatic corrections in the path 
integral formulation can be found in Ref. 101 . For each r the influence functional 
Z[R] is given by the partition function of the electronic subsystem evaluated for 
the respective nuclear configuration R(r) . Assuming that the temperature is small 
compared to the gap in the electronic spectrum only the electronic ground state 
with energy E (R(t)) (obtained from solving Eq. (20) without the internuclear 
Coulomb repulsion term) is populated. This electronic ground state dominance 
leads to the following simple expression 



Z[R]bo = cxp 



dTE (R(T)) 



which yields the final result 



Z BO - j^VRcxp - dr (T(R) + V(K) + £ (R)) 



(313) 



(314) 



Here nuclear exchange is neglected by assuming that the nuclei are distinguishable 
so that they can be treated within Boltzmann statistics, which corresponds to the 
Hartree approximation for the nuclear density matrix. The presentation given here 
follows Ref. 399 and alternative derivations were given in Sect. 2.3 of Refs. 124 and 
in the appendix of Ref. 427 . There, a wavefunction basis instead of the position 
basis as in Eq. (312) was formally used in order to evaluate the influence functional 
due to the electrons. 

The partition function Eq. (314) together with the Coulomb Hamiltonian Eq. (2) 
leads after applying the lowest-order Trotter factorization 334 to the following dis- 
cretized expression 

p N 

Zb ° = n n 



s=l 1=1 



s) 



x exp 



-0 



t if \ m ^p K - *4 S+1) ) 2 + > ({R,} (s) ) 



s=i yi=i 

,2 _ -piai 



(315) 



for the path integral with tup = P/0 2 Thus, the continuous parameter r S [0,0] is 
discretized using P so-called Trotter slices or "time slices" s = 1 , . . . , P of "dura- 
tion" At = 0/P. The paths 



{{R 7 } (S) } = ({R/} (1) ;...;{R/} (P) ) 



R (l). .T>(P) 



,R 



(P)> 



(316) 



have to be closed due to the trace condition, i.e. they are periodic in imaginary 
time r which implies R/(0) = R/(/3) and thus Rj 



r', 1 '; the internuclear 



Coulomb repulsion T^(R) is now included in the definition of the total electronic 
energy E - Note that Eq. (315) is an exact reformulation of Eq. (314) in the limit 
of an infinitely fine discretization P — > oo of the paths. 
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The effective classical partition function Eq. (315) with a fixed discretization P 
is isomorphic to that for N polymers each comprised by P monomers 23 3,i26,i20_ 
Each quantum degree of freedom is found to be represented by a ring polymer or 
necklace. The intrapolymeric interactions stem from the kinetic energy T(R) and 
consist of harmonic nearest neighbor couplings otwp along the closed chain. The 
interpolymeric interaction is given by the scaled potential / P which is only 
evaluated for configurations {R/}^ at the same imaginary time slice s. 

In order to evaluate operators based on an expression like Eq. (315) most nu- 
merical path integral schemes utilize Metropolis Monte Carlo sampling with the 
effective potential 



\ Ml ul (r« - R< s+1) ) 2 + p£o ({R/} (s) ) } (317) 





of the isomorphic classical system 233,i26,542,i20,i24,646,407_ Molecular dynam- 
ics techniques were also proposed in order to sample configuration space, see 
Refs. 99 > 490 > 4 62, 501,273 f or pi one ering work and Ref. 646 for an authoritative review. 
Formally a Lagrangian can be obtained from the expression Eq. (317) by extending 
it 

" 't - 5"** K - R i' +1 ') 2 ) - ({*<>"> 

1=1 v 1 7 

(318) 

with NxP fictitious momenta P ^ and corresponding (unphysical) fictitious masses 
M'j. At this stage the time dependence of positions and momenta and thus the time 
evolution in phase space as generated by Eq. (318) has no physical meaning. The 
sole use of "time" is to parameterize the deterministic dynamical exploration of 
configuration space. The trajectories of the positions in configuration space, can, 
however, be analyzed similar to the ones obtained from the stochastic dynamics 
that underlies the Monte Carlo method. 

The crucial ingredient in ah initio 395,399,644,404 as opposed to stan- 
dard 233 . !26, 542, 120, 124, 646, 407 p a th integral simulations consists in computing the 
interactions Eq "on-the-fly" like in ah initio molecular dynamics. In analogy to this 
case both the Car-Parrinello and Born-Oppcnhcimcr approaches from Sects. 2.4 
and 2.3, respectively, can be combined with any electronic structure method. The 
first implementation 395 was based on the Car-Parrinello / density functional com- 
bination from Sect. 2.4 which leads to the following extended Lagrangian 



,{R/} 




R< s+1) ) ) , (319) 
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where the interaction energy E KS [{fa}^, {Rj}W] at time slice s is denned in 
Eq. (75); note that here and in the following the dots denote derivatives with respect 
to propagation time t and that E^ s = min £; KS . The standard Car-P arrincllo 
Lagrangian, see e.g. Eq. (41) or Eq. (58), is recovered in the limit P — 1 which 
corresponds to classical nuclei. Mixed classical / quantum systems can easily be 
treated by representing an arbitrary subset of the nuclei in Eq. (319) with only one 
time slice. 

This simplest formulation of ab initio path integrals, however, is insufficient for 
the following reason: ergodicity of the trajectories and adiabaticity in the sense 
of Car-Parrinello simulations are not guaranteed. It is known since the very first 
molecular dynamics computer experiments that quasiharmonic systems (such as 
coupled stiff harmonic oscillators subject to weak anharmonicities, i.e. the famous 
Fermi-Pasta-Ulam chains) can easily lead to nonergodic behavior in the sampling 
of phase space 210 . Similarly "microcanonical" path integral molecular dynamics 
simulations might lead to an insufficient exploration of configuration space depend- 
ing on the parameters 273 . The severity of this nonergodicity problem is governed 
by the stiffness of the harmonic intrachain coupling oc ujp and the anharmonicity of 
the overall potential surface oc E KS /P which establishes the coupling of the modes. 
For a better and better discretization P the harmonic energy term dominates ac- 
cording to ~ P whereas the mode-mixing coupling decreases like ~ 1/P. This 
problem can be cured by attaching Nose-Hoover chain thermostats 388 , see also 
Sect. 4.2, to all path integral degrees of freedom 637 > 644 . 

The second issue is related to the separation of the power spectra associated 
to nuclear and electronic subsystems during Car-Parrinello ab initio molecular dy- 
namics which is instrumental for maintaining adiabaticity, see Sect. 2.4. In ab 
initio molecular dynamics with classical nuclei the highest phonon or vibrational 
frequency a;™ 3 ^ is dictated by the physics of the system, see e.g. Fig. 2. This means 
in particular that an upper limit is given by stiff intramolecular vibrations which 
do not exceed w™ ax < 5000 cm -1 or 150 THz. In ab initio path integral simula- 
tions, on the contrary, w™ ax is given by top which actually diverges with increasing 
discretization as ~ \fP. The simplest counteraction would be to compensate this 
artifact by decreasing the fictitious electron mass \i until the power spectra are 
again separated for a fixed value of P and thus ujp. This, however, would lead to 
a prohibitively small time step because Ai max oc ^JJL. This dilemma can be solved 
by thermostatting the electronic degrees of freedom as well 395 > 399 > 644 ; se e Sect. 4.2 
for a related discussion in the context of metals. 

Finally, it is known that diagonalizing the harmonic spring interaction in 
Eq. (319) leads to more efficient propagators 637 > 644 . One of these transforma- 
tion and the resulting Nose-Hoover chain thermostatted equations of motion will 
be outlined in the following section, see in particular Eqs. (331)-(337). In addi- 
tion to keeping the average temperature fixed it is also possible to generate path 
trajectories in the isobaric-isothermal NpT ensemble 646 > 392 . Instead of using 
Car-Parrinello fictitious dynamics in order to evaluate the interaction energy in 
Eq. (318), which is implemented in the CPMD package 142 , it is evident that also 
the Born-Oppenheimer approach from Sect. 2.3 or the free energy functional from 
Sect. 4.3 can be used. This route eliminates the adiabaticity problem and was taken 
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Up e.g. in Rcfs. 132,37,596,597,428,429,333 



A final observation concerning parallel supercomputers might be useful, see also 
Sect. 3.9. It is evident from the Lagrangian Eq. (319) and the resulting equations 
of motion (e.g. Eqs. (331)-(337)) that most of the numerical workload comes from 
calculating the ab initio forces on the nuclei. Given a fixed path configuration 
Eq. (316) the P underlying electronic structure problems are independent from 
each other and can be solved without communication on P nodes of a distributed 
memory machine. Communication is only necessary to send the final result, essen- 
tially the forces, to a special node that computes the quantum kinetic contribution 
to the energy and integrates finally the equations of motions. It is even conceiv- 
able to distribute this task on different supercomputers, i.e. "meta-computing" is 
within reach for such calculations. Thus, the algorithm is "embarrassingly parallel" 
provided that the memory per node is sufficient to solve the complete Kohn-Sham 
problem at a given time slice. If this is not the case the electronic structure cal- 
culation itself has to be parallelized on another hierarchical level as outlined in 
Sect. 3.9. 

4-4-3 Ab Initio Path Centroids: Dynamics 

Initially the molecular dynamics approach to path integral simulations was invented 
merely as a trick in order to sample configuration space similar to the Monte Carlo 
method. This perception changed recently with the introduction of the so-called 



"centroid molecular dynamics" technique 102 , see Refs. 103,104,105,665,505,506,507 f or 



background information. In a nutshell it is found that the time evolution of the 
centers of mass or centroids 



of the closed Feynman paths that represent the quantum nuclei contains quasi- 
classical information about the true quantum dynamics. The centroid molecular 
dynamics approach can be shown to be exact for harmonic potentials and to have 
the correct classical limit. The path centroids move in an effective potential which 
is generated by all the other modes of the paths at the given temperature. This 
effective potential thus includes the effects of quantum fluctuations on the (qua- 
siclassical) time evolution of the centroid degrees of freedom. Roughly speaking 
the trajectory of the path centroids can be regarded as a classical trajectory of the 
system, which is approximately "renormalized" due to quantum effects. 

The original centroid molecular dynamics technique 102 , 103,104,105, 665 re ]j es on 
the use of model potentials as the standard time independent path integral simu- 
lations. This limitation was overcome independently in Refs. 469 ' 411 by combining 
ab initio path integrals with centroid molecular dynamics. The resulting technique, 
ab initio centroid molecular dynamics can be considered as a quasiclassical gener- 
alization of standard ab initio molecular dynamics. At the same time, it preserves 
the virtues of the ab initio path integral technique 395 > 3 "> 644 ' 404 to generate exact 
time-independent quantum equilibrium averages. 




(320) 



s 



'=1 
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Here, the so-called adiabatic formulation 105 ' 390 > 106 f a f) i n Mo centroid molecu- 
lar dynamics 411 is discussed. In close analogy to ab initio molecular dynamics with 
classical nuclei also the effective centroid potential is generated "on-the-fly" as the 
centroids are propagated. This is achieved by singling out the centroid coordinates 
in terms of a normal mode transformation 138 and accelerating the dynamics of 
all non-centroid modes artificially by assigning appropriate fictitious masses. At 
the same time, the fictitious electron dynamics a la Car-Parrinello is kept in order 
to calculate efficiently the ab initio forces on all modes from the electronic struc- 
ture. This makes it necessary to maintain two levels of adiabaticity in the course 
of simulations, see Sect. 2.1 of Ref. 411 for a theoretical analysis of that issue. 

The partition function Eq. (315), formulated in the so-called "primitive" path 
variables {Rj}( s ), is first transformed 644 > 646 to a representation in terms of the 
normal modes {u/}^\ which diagonalizc the harmonic nearest-neighbor harmonic 
coupling 138 . The transformation follows from the Fourier expansion of a cyclic 
path 

p 

= a/" ex P i 2m ( s - !)( s ' - l )l p \ > ( 321 ) 

s'=l 

where the coefficients {a/}( s ) are complex. The normal mode variables {uj}( s ) are 
then given in terms of the expansion coefficients according to 

„(P) _ _((^+2)/2) 

uf- 2) =Re (a< s) ) 

uf s - 1] =Im (a< s) ) . (322) 

Associated with the normal mode transformation is a set of normal mode frequencies 
{A}( s ) given by 



^(28-1) = ^(28-2) = 2 p 



1 — cos 



/2tt(s - 1) 

V P 



(323) 



with A^ = and A( p ) = 4P. Equation (321) is equivalent to direct diagonalization 
of the matrix 

A ss / = 2S SS ' - <5 s ,s'-i - ^s,s'+i (324) 

with the path periodicity condition A s o = A s p and A S; p_|_i = A s i and subsequent 
use of the unitary transformation matrix U to transform from the "primitive" 
variables {R/}( s ) to the normal mode variables {uj}^) 

R^.vpfuLur) 

s' = l 

p 
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The eigenvalues of A when multiplied by P are precisely the normal mode frequen- 
cies {A}( s ). Since the transformation is unitary its Jacobian is unity. Finally it is 
convenient to define a set of normal mode masses 

M\ s) = A (s) M/ (326) 

that vary along the imaginary time axis s — 1,...,P, where A^ 1 ^ = for the 
centroid mode . 

Based on these transformations the Lagrangian corresponding to the ab initio 
path integral expressed in normal modes is obtained 644 



ipi 



{*} W ,{R/ (u« ...,ur)} 




(327) 

\ / \ / I 

1=1 ) 

where the masses M'j 8 ^ will be defined later, see Eq. (338). As indicated, the 
electronic energy E^ is always evaluated in practice in terms of the "primitive" 
path variables {R/}' s ^ in Cartesian space. The necessary transformation to switch 
forth and back between "primitive" and normal mode variables is easily performed 
as given by the relations Eq. (325). 

The chief advantage of the normal mode representation Eq. (325) for the present 
purpose is that the lowest-order normal mode u^ 1 ' 

u« = = i f; R<''> (328) 

s'=l 

turns out to be identical to the centroid Rj of the path that represents the 7th 
nucleus. The centroid force can also be obtained from the matrix U according 



to 644 



dE 1 .A dE (s "> 



%r " " (329) 



since Ui a = Ujj = 1/VF and the remaining normal mode forces are given by 

—— p- = —j= U ss / — — r-pr for« = 2,...,P (330) 

in terms of the "primitive" forces —dE^/dH^. Here, E on the left-hand-side 
with no superscript (s) refers to the average electronic energy E = (1 / P) J2s=i 
from which the forces have to be derived. Thus, the force Eq. (329) acting on each 
centroid variable u^ 1 ', I = 1,...N, is exactly the force averaged over imaginary 
time s — 1, . . . , P, i.e. the centroid force on the Ith nucleus as already given in 
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Eq. (2.21) of Ref. 644 . This is the desired relation which allows in centroid molecular 
dynamics the centroid forces to be simply obtained as the average force which acts 
on the lowest-order normal mode Eq. (328). The non-centroid normal modes 
s = 2, 3, . . . , P of the paths establish the effective potential in which the centroid 
moves. 

At this stage the equations of motion for adiabatic ab initio centroid molecular 
dynamics 411 can be obtained from the Euler-Lagrange equations. These equations 
of motion read 



= 4£ 
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where denotes the Cartesian components of a given normal mode vector 

vSf 1 = (uj S \, uf\, u\ s \). In the present scheme, independent Nose-Hoover chain 
thermostats 388 of length K are coupled to all non-centroid mode degrees of free- 
dom s = 2, . . . , P 
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and all orbitals at a given imaginary time slice s are thermostatted by one such 
thermostat chain of length L 



Qlvi s) = 2 



Em 
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note that for standard ab initio path integral runs as discussed in the previous 
section the centroid mode should be thermostatted as well. The desired fictitious 
kinetic energy of the electronic subsystem T® can be determined based on a short 
equivalent classical Car-Parrinello run with P = 1 and using again the relation 
= 2T®/QN' e where N' e is the number of orbitals. The mass parameters {Q1} 
associated to the orbital thermostats are the same as those defined in Eq. (271), 
whereas the single mass parameter Q n for the nuclei is determined by the harmonic 
interaction and is given by Q n = k^T/tup = /3/P. The characteristic thermostat 



116 



frequency of the electronic degrees of freedom uj c should again lie above the fre- 
quency spectrum associated to the fictitious nuclear dynamics. These is the method 
that is implemented in the CPMD package 142 . 

An important issue for adiabatic ab initio centroid molecular dynamics 411 is 
how to establish the time-scale separation of the non-centroid modes compared to 
the centroid modes. This is guaranteed if the fictitious normal mode masses 
are taken to be 

Mj (1) = Mi 

M; (s) = 7 Mj s) ,s = 2,...,P , (338) 

where Mj is the physical nuclear mass, are the normal mode masses Eq. (326), 
and 7 is the "centroid adiabaticity parameter" ; note that this corrects a misprint 
of the definition of Mj for s > 2 in Ref. 411 . By choosing 0(7 <C 1, the required 
time-scale separation between the centroid and non-centroid modes can be con- 
trolled so that the motion of the non-centroid modes is artificially accelerated, see 
Sect. 3 in Ref. 411 for a systematic study of the 7-dependence. Thus, the centroids 
with associated physical masses move quasiclassically in real time in the centroid 
effective potential, whereas the fast dynamics of all other nuclear modes s)l is fic- 
titious and serves only to generate the centroid effective potential "on-the-fly" . In 
this sense 7 (or rather 7M/) is similar to fi, the electronic adiabaticity parameter 
in Car-Parrinello molecular dynamics. 



4-4-4 Other Approaches 

It is evident from the outset that the Born-Oppcnheimcr approach to generate 
the ab initio forces can be used as well as Car-Parrinello molecular dynamics 
in order to generate the ab initio forces on the quantum nuclei. This varia- 
tion was utilized in a variety of investigations ranging from clusters to molecular 
solids !32, 37,596, 597,428, 429,333^ closely related to the ab initio path integral approach 
as discussed here is a method that is based on Monte Carlo sampling of the path 
integral 672 . It is similar in spirit and in its implementation to Born-Oppcnheimcr 
molecular dynamics sampling as long as only time-averaged static observables are 
calculated. A scmicmpirical ("cndo" and "indo" ) version of Born-Oppcnhcimer ab 
initio path integral simulations was also devised 656 and applied to study muonated 
organic molecules 656 > 657 . 

A non-sclf-consistent approach to ab initio path integral calculations was advo- 
cated and used in a series of publications devoted to study the interplay of nuclear 
quantum effects and electronic structure in unsaturated hydrocarbons like ben- 
zene 544 ! 503 ,8i,543,504^ According to this philosophy, an ensemble of nuclear path 
configurations Eq. (316) is first generated at finite temperature with the aid of a 
parameterized model potential (or using a tight-binding Hamiltonian 504 ). In a sec- 
ond, independent step electronic structure calculations (using Pariser-Parr-Poplc, 
Hubbard, or Hartrcc-Fock Hamiltonians) arc performed for this fixed ensemble of 
discretized quantum paths. The crucial difference compared to the self consistent 
approaches presented above is that the creation of the thermal ensemble and the 
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subsequent analysis of its electronic properties is performed using different Hamil- 
tonians. 

Several attempts to treat also the electrons in the path integral formulation - 
instead of using wavefunctions as in the ab initio path integral family - were 
published 606,119,488, 449,450^ These approaches are exact in principle, i.e. non 
adiabaticity and full electron phonon coupling is included at finite temperatures. 
However, they suffer from severe stability problems 121 in the limit of degenerate 
electrons, i.e. at very low temperatures compared to the Fermi temperature, which 
is the temperature range of interest for typical problems in chemistry and materials 
science. Recent progress on computing electronic forces from path integral Monte 
Carlo simulations was also achieved 708 . 

More traditional approaches use a wavefunction representation for both the 
electrons in the ground state and for nuclear density matrix instead of path in- 
tegrals. The advantage is that real-time evolution is obtained more naturally 
compared to path integral simulations. A review of such methods with the em- 
phasis of computing the interactions "on-the-fly" is provided in Ref. 158 . An ap- 
proximate wavefunction-based quantum dynamics method which includes several 
excited states and their couplings was also devised and used 385 > 386 >387,45^ ^ n a ^ er _ 
native approach to approximate quantum dynamics consists in performing instan- 
ton or semiclassical ab initio dynamics 325 ' 47 . Also the approximate vibrational 
self-consistent field approach to nuclear quantum dynamics was combined with 
"on-the-fly" MP2 electronic structure calculations 122 . 



5 Applications: From Materials Science to Biochemistry 



5. 1 Introduction 

Ab initio molecular dynamics was called a "virtual matter laboratory" 234 , a notion 
that is fully justified in view of its relationship to experiments performed in the 
real laboratory. Ideally, a system is prepared in some initial state and than evolves 
according to the basic laws of physics - without the need of experimental input. 
It is clear to every practitioner that this viewpoint is highly idealistic for more 
than one reason, but still this philosophy allows one to compute observables with 
predictive power and also implies a broad range of applicability. 

It is evident from the number of papers dealing with ab initio molecular dy- 
namics, see for instance Fig. 1, that a truly comprehensive survey of applications 
cannot be given. Instead, the strategy chosen is to provide the reader with a wealth 
of references that try cover the full scope of this approach - instead of discussing 
in depth the physics or chemistry of only a few specific applications. To this end 
the selection is based on a general literature search in order to suppress personal 
preferences as much as possible. In addition the emphasis lies on recent applica- 
tions that could not be covered in earlier reviews. This implies that several older 
key reference papers on similar topics are in general missing. 
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5.2 Solids, Polymers, and Materials 



The first application of Car-Parrinello molecular dynamics 108 dealt with silicon, 
one of the basic materials in semiconductor industry. Classic solid-state applica- 
tion of this technique focus on the properties of crystals, such as those of CuCl 
where anharmonicity and off-center displacements of the Cu along the (111) di- 
rections were found to be important to describe the crystal structure as a func- 
tion of temperature and pressure 64 . Various properties of solid nitromethane 647 , 
crystalline nitric acid trihydrate 602 , solid benzene 420 , stage-1 alkali-graphite in- 
tercalation compounds 286 > 287 ; and of the one-dimensional intercalation compound 
2HgS»SnBr2 530 were determined based on first principles. The molecular solid HBr 
undergoes various phase transitions upon compression. The dynamical behavior of 
one of these phases, disordered HBr-I, could be clarified using ab initio molecular 
dynamics 313 . Structure, phase transitions and short time dynamics of magnesium 
silicate perovskites were analyzed in terms of ab initio trajectories 670 . The A7 to 
simple cubic transformation in As was investigated using ab initio molecular dynam- 
ics at constant-pressure 568 . By applying external pressure the hydrogen sublattice 
was found to undergo amorphization in Mg(OH) 2 and Ca(OH) 2 a phenomenon 
that was interpreted in terms of frustration 511 . Properties of solid cubane CsHs 
were obtained in constant pressure simulations and compared to experiment 514 . 
Ab initio simulations of the graphitization of flat and stepped diamond (111) sur- 
faces uncovered that the transition temperature depends sensibly on the type of 
the surface 327 . 

Sliding of grain boundaries in aluminum as a typical ductile metal was generated 
and analyzed in terms of atomistic rearrangements 432 . Microfracture in a sample 
of amorphous silicon carbide was induced by uniaxial strain and found to induce Si 
segregation at the surface 226 . The early stages of nitride growth on cubic silicon 
carbide including wetting were modeled by depositing nitrogen atoms on the Si- 
terminated SiC(OOl) surface 225 . 

Classical proton diffusion in crystalline silicon at high temperatures was an early 
application to the dynamics of atoms in solids 93 . Using the ab initio path integral 
technique 395 . 3 "i 644 » 404 ) se e Sect. 4.4 the preferred sites of hydrogen and muonium 
impurities in crystalline silicon 428 > 429 ; or the proton positions in HC1 • nH 2 crys- 
talline hydrates 516 could be located. The radiation-induced formation of defects 
in c-Si via vacancies and self intcrstitials was simulated by ab initio molecular dy- 
namics 178 . The classical diffusion of hydrogen in crystalline GaAs was followed in 
terms of diffusion paths 668 and photoassisted reactivation of H-passivated Si donors 
in GaAs was simulated based on first principles 430 . Oxygen diffusion in p-doped 
silicon can be enhanced by adding hydrogen to the material, an effect that could 
be rationalized by simulations 107 . Ab initio dynamics helped in quantifying the 
barrier for the diffusion of atomic oxygen in a model silica host 279 . The microscopic 
mechanism of the proton diffusion in protonic conductors, in particular Sc-doped 
SrTiOa and Y-doped SrCe03, is studied via ab initio molecular dynamics, where 
is it found that covalent OH-bonds are formed during the process 561 . Ionic dif- 
fusion in a ternary superionic conductor was obtained by ab initio dynamics 677 . 
Proton motion and isomerization pathways of a complex photochromic molecular 
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crystal composed of 2-(2,4-dinitrobenzyl)pyridine dyes was generated by ab initio 
methods 216 . 

Also materials properties of polymers are investigated in quite some detail. 
Early applications of scmicmpirical zdo molecular dynamics 666 were devoted to 
defects in conducting polymers, in particular to solitons, polarons and alkali doping 
in polyacetylene 666 ' 667 as well as to muonium implanted in trans and cis polyacety- 
lene 200 . More recent are calculations of Young's modulus for crystalline polyethy- 
lene 271 , soliton dynamics in positively charged polyacetylene chains 125 , charge 
localization in doped polypyrroles 140 , chain rupture of polyethylene chains under 
tensile load 533 , the influence of a knot on the strength of a polymer strand 534 , or 
ion diffusion in polyethylene oxide 456 . 



5.3 Surfaces, Interfaces, and Adsorbates 

A host of studies focusing on atoms and in particular on molecules interacting 
with surfaces appeared over the years. Recent studies focussed for instance on 
C2H2, C2H4, and trimethylgallium adsorbates on the GaAs(001)-(2x4) surface 248 , 
thiophcnc on the catalytically active MoS2(010) 512 or RuS2 580 surfaces, small 
molecules on a nitric acid monohydrate crystal surface 624 , CO on Si(OOl) 314 , small 
molecules on Ti0 2 554 ' 41 , sulfur on Si(100) at various coverages , and sulfuric 
acid adsorbed on ZrO 2 (101) and ZrO 2 (001) 269 . 

Specific to ab initio molecular dynamics is its capability to describe also 
chcmisorption as well as dynamical processes on (and of) surfaces including surface 
reactions 500 . The ab initio calculations of surface phonons in semiconductor sur- 
faces can be based on the frozcn-phonon, linear-response or nowadays molecular 
dynamics approaches, see Ref. 218 for a discussion and comparison. A review on 
the structure and energetics of oxide surfaces including molecular processes occur- 
ring on such surfaces is provided in Ref. 235 , whereas Ref. 256 concentrates on the 
interaction of hydrogen with clean and adsorbate covered metal and semiconductor 
surfaces. 

Recent applications in surface science include the transition from surface vibra- 
tions to liquid-like diffusional dynamics of the Ge(lll) surface 607 , the diffusion of Si 
adatoms on a double-layer stepped Si(OOl) surface 33 °, the structure of chemisorbed 
acetylene on the Si(001)-(2x 1) surface 423 , chemisorption of quinizarin on a 
AI2O3 212 > 213 ; the diffusion of a single Ga adatom on the GaAs(100)-c(4x4) sur- 
face 367 , homoepitaxial crystal growth on Si(OOl) and the low-temperature dynam- 
ics of Si(lll)-(7x7) 595 ' 611 , dissociation of an H 2 molecule on MgO 358 > 359 , disso- 
ciation of CI2 on GaAs(llO) 3S0 , chlorine adsorption and reactions on Si(100) 691 , 
molecular motion of NH 3 on MgO 358 , dynamics and reactions of hydrated a- 
alumina surfaces 289 , molecular vs. dissociative adsorption of water layers on 
MgO(lOO) as a function of coverage 448 , oxidation of CO on Pt(lll) 8 < 705 , the 
reaction HC1 + HOC1 — > H 2 + CI2 as it occurs on an ice surface 373 , or desorp- 
tion of D2 from Si(100) 255 . Thermal contraction, the formation of adatom-vacancy 
pairs, and finally premelting was observed in ab initio simulations of the Al(llO) 
surface at temperatures up to 900 K 415 Early stages of the oxidation of a Mg(0001) 
surface by direct attack of molecular O2 was dynamically simulated 96 including 
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the penetration of the oxidation layer into the bulk. Similarly, the growth of an 
oxide layer was generated on an Si(100) surface 653 . 

The water-Pd(lOO), water-O/Pd(100) and water-Si(lll) interfaces were simu- 
lated based on ab initio molecular dynamics 336 > 655 . Water covering the surface of 
a microscopic model of muscovite mica is found to form a two-dimensional network 
of hydrogen bonds, called two-dimensional ice, on that surface 447 . The metal- 
organic junction of monolayers of Pd-porphyrin and perylene on Au(lll) was ana- 
lyzed using an ab initio approach 355 . An interesting possibility is to compute the 
tip-surface interactions in atomic force microscopy as e.g. done for a neutral silicon 
tip interacting with an InP(llO) surface 619 or Si(lll) 481 < 482 . 



5.4 Liquids and Solutions 

Molecular liquids certainly belong to the classic realm of molecular dynamics simu- 
lations. Water was and still is a challenge 581 for both experiment and simulations 
due to the directional nature and the weakness of the hydrogen bonds which leads 
to delicate association phenomena. Pioneering ab initio simulations of pure wa- 
ter at ambient 352 and supercritical conditions 205 were reported only a few years 
ago. More recently, these gradient-corrected density functional theory-based simu- 
lations were extended into several directions 587 > 573 < 575 > 576 ,579,ii8^ j n ^ e mean time 
(minimal-basis) Hartree-Fock ab initio molecular dynamics 291 as well as more ap- 
proximate schemes 455 were also applied to liquid water. Since chemical reactions 
often occur in aqueous media the solvation properties of water are of utmost impor- 
tance so that the hydration of ions 403,620,621,377,502 and gmaU mo i ccmes 353,354,433 

was investigated. Similarly to water liquid HF is a strongly associated liquid which 
features short-lived hydrogen-bonded zig-zag chains 521 . Another associated liq- 
uid, methanol, was simulated at 300 K using an adaptive finite-element method 634 
in conjunction with Born-Oppcnheimer molecular dynamics 635 . In agreement with 
experimental evidence, the majority of the molecules is found to be engaged in short 
linear hydrogen-bonded chains with some branching points 635 . Partial reviews on 
the subject of ab initio simulations as applied to hydrogen-bonded liquids can be 
found in the literature 586,406,247^ 

The ab initio simulated solvation behavior of "unbound electrons" in liquid 
ammonia at 260 K was found to be consistent with the physical picture extracted 
from experiment 155 > 156 . Similarly, ab initio molecular dynamics of dilute 553 > 203 and 
concentrated 569 molten K x -(KCl)i_ x mixtures were performed at 1300 K entering 
the metallic regime. The structure of liquid ammonia at 273 K was investigated 
with a variety of techniques so that limitations of using classical nuclei, simple point 
charge models, small systems, and various density functional could be assessed 164 . 

Ab initio molecular dynamics is also an ideal tool to study other complex fluids 
with partial covalcncy, metallic fluids, and their transformations as a function of 
temperature, pressure, or concentration. The properties of water-free KF • nHF 
melts depend crucially on polyfluoride anions H m F~ +1 and solvated K + cations. 
Ab initio simulations allow for a direct comparison of these complexes in the liquid, 
gaseous and crystalline phase 515 . The changes of the measured structure factor 
of liquid sulfur as a function of temperature can be rationalized on the atomistic 
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level by various chain and ring structures that can be statistically analyzed in ab 
initio molecular dynamics simulations 631 . Liquid GeSe2 is characterized by strong 
chemical bonds that impose a structure beyond the usual very short distances due 
to network formation 416 . Zintl-alloys such as liquid NaSn 552 or KPb 556 have very 
interesting bonding properties that manifest themselves in strong temperature- and 
concentration dependences of their structure factors (including the appearance of 
the so-called first sharp diffraction peak 555 ) or electric conductivities. 

Metals are ideal systems to investigate the metal-insulator transition upon ex- 
pansion of the liquid 346 ' 63 or melting 689 . Liquid copper was simulated at 1500 K: 
structural and dynamical data were found to be in excellent agreement with exper- 
imental 464 . Transport coefficients of liquid metals (including in particular extreme 
conditions) can also be obtained from first principles molecular dynamics using the 
Green-Kubo formalism 571 ' 592 . The microscopic mechanism of the semiconductor- 
metal transition in liquid As2Se 3 could be rationalized in terms of a structural 
change as found in ab initio simulations performed as a function of temperature 
and pressure 563 . The III— V semiconductors, such as GaAs, assume metallic behav- 
ior when melted, whereas the II- VI semiconductor CdTe does not. The different 
conductivities could be traced back to pronounced structural dissimilarities of the 
two systems in the melt 236 . 



5.5 Glasses and Amorphous Systems 

Related to the simulation of dynamically disordered fluid systems are investiga- 
tions of amorphous or glassy materials. In view of the severe limitations on system 
size and time scale (and thus on correlation lengths and times) ab initio molecular 
dynamics can only provide fairly local information in this sense. Within these inher- 
ent constraints the microscopic structure of amorphous selenium 304 and tetrahedral 
amorphous carbon 384 , the amorphization of silica 684 , boron doping in amorphous 
Si:H 181 or in tetrahedral amorphous carbon 227 , as well as the Raman spectrum 465 
and dynamic structure factor 466 of quartz glass and their relation to short-range 
order could be studied. 

The properties of supercooled CdTe were compared to the behavior in the liq- 
uid state in terms of its local structure 237 . Defects in amorphous Sii^Ge^ alloys 
generated by ab initio annealing were found to explain ESR spectra of this sys- 
tem 329 . The infrared spectrum of a sample of amorphous silicon was obtained 
and found to be in quantitative agreement with experimental data 152 . The CO2 
insertion into a model of argon-bombarded porous Si02 was studied 508 . In partic- 
ular the electronic properties of amorphous GaN were investigated using ab initio 
methods 601 . 

Larger systems and longer annealing times are accessible after introducing more 
approximations into the first principle treatment of the electronic structure that 
underlies ab initio molecular dynamics. Using such methods 551 , a host of different 
amorphous carbon nitride samples with various stoichiometrics and densities could 
be generated and characterized in terms of trends 675 . Similarly, the pressure- 
induced glass-to-crystal transition in condensed sodium was investigated 22 and two 
structural models of amorphous GaN obtained at different densities were examined 
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in terms of their electronic structure 



5.6 Matter at Extreme Conditions 

A strong advantage of ab initio simulations is their predictive power also at ex- 
treme conditions, an area where molecular dynamics relying on fitted potential 
models might encounter severe difficulties. Thus, high pressures and / or high tem- 
peratures such as those in the earth's core, on other planets, or on stars can be 
easily achieved in the virtual laboratory. This opens up the possibility to study 
phase transformations and chemical reactions at such conditions 56 . Furthermore, 
conditions of geophysical and astrophysical interest can nowadays be produced in 
the real laboratory, using techniques based on diamond anvil cells, shock waves, or 
lasers. The limitations of these experimental approaches are, however, not so much 
related to generating the extreme conditions as one might expect, but rather to 
measuring observables. 

In the virtual laboratory this information is accessible and the melting of 
diamond at high pressure 222 , the phase transformation from the antiferromag- 
netic insulating 5-02 phase to a nonmagnetic metallic molecular C~C>2 phase 557 , 
the phase diagram of carbon at high pressures and temperatures 261 as well as 
transformations of methane 13 , carbon monoxide 54 , molecular CO2 267 > 558 , water 

ice 363,364,58,50,51,52 5 soM 305,337,65,66,333 and hot fluid 5 hydrogcn; solid Ar(H 2 ) 2 53 

under pressure could be probed. Along similar lines properties of a liquid Fe-S mix- 
ture under earth's core conditions 11 , the viscosity of liquid iron 690 > 592 ; the sound 
velocity of dense hydrogen at conditions on jupiter 6 , the phase diagram of water 
and ammonia up to 7000 K and 300 GPa 118 , the laser heating of silicon 570 > 572 
and graphite 574 etc. were investigated at extreme state points. A review on ab 
initio simulations relevant to minerals at conditions found in the earth's mantle is 
provided in Rcf . 683 . 

5.7 Clusters, Fullerenes, and Nanotubes 

Investigations of clusters by ab initio molecular dynamics were among the first 
applications of this technique. Here, the feasibility to conduct finite-temperature 
simulations and in particular the possibility to search globally for minima turned 
out to be instrumental 302 . 31 . 303 > 550 » 517 . 519 ) S ee e.g. Refs. 16 > 321 > 32 for reviews. 
Such investigations focus more and more on clusters with varying composi- 
tion s 18 ' 293 . 199 ' 348 ' 349 * 16 ^ Cluster melting is also accessible on an ab initio foot- 
ing 84,531,525,526 and mo i ecmar clusters, complexes or cluster aggregates are actively 
investigated 612,645, 613, 70,596, 597,133, 701,524 

In-v semiconductor clusters embedded in sodalite show quantum confinement 
and size effects that can be rationalized by ab initio simulations 625 > 95 . Supported 
clusters such as Cu„ on an MgO(100) surface are found to diffuse by "rolling" and 
"twisting" motions with very small barriers 438 . The diffusion of protonated helium 
clusters in various sodalite cages was generated using ab initio dynamics 198 . Photo- 
induced structural changes in Se chains and rings were generated by a vertical 
homo — * lumo excitation and monitored by ab initio dynamics 306 . With the 
discovery and production of finite carbon assemblies ab initio investigations of the 
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properties of fullcrenes 19 > 17 > 451 ; the growth process of nanotubes 127 > 62 > 72 ; or the 
electrical conductivity of nanowires 38 > 272 became of great interest. 

5.8 Complex and Floppy Molecules 

The determination of the structure of a RNA duplex including its hydration wa- 
ter 311 , investigations of geometry and electronic structure of porphyrins and por- 
phyrazines 356 , and the simulation of a bacteriochlorophyll crystal 381 are some 
applications to large molecules. Similarly the "carboplatin" complex 623 - a drug 
with large ligands - as well as the organometallic complex Alq(3) 148 - an electro- 
luminescent material used in organic light emitting diodes - were investigated with 
respect to structural, dynamical and electronic properties. 

The organometallic compound C2H2Li2 has an unexpected ground-state struc- 
ture that was found by careful ab initio simulated annealing 521 . In addition, this 
complex shows at high temperatures intramolecular hydrogen migration that is 
mediated via a lithium hydride unit 521 . Ground-state fluxionality of protonated 
methane CH^" 397 ' 408 including some isotopomers 409 and of protonated acetylene 
C2H3" 400 was shown to be driven by quantum effects. The related dynamical 
exchange of atoms in these molecules can also be excited by thermal fluctua- 
tions 630 > 85 > 401 . in addition it was shown that CH^" is three-center two-electron 
bonded and that this bonding topology does not qualitatively change in the pres- 
ence of strong quantum motion 402 . The fluxional behavior of the protonated ethane 
molecular ion C2H7 was investigated by ab initio molecular dynamics as well 172 . 

The neutral and ionized SiHs and Si2H3 species display a rich dynamical be- 
havior which was seen during ab initio molecular dynamics simulations 246 . The 
lithium pentamer Lis wa s found to perform pseudorotational motion on a time 
scale of picoseconds or faster at temperatures as low as 77 K 231 . Using ab initio 
instanton dynamics the inversion splitting of the NH 3 , ND 3 , and PH 3 molecules 
due to the umbrella mode was estimated 325 . Similarly, a semiclassical ab initio 
dynamics approach as used to compute the tunneling rate for intramolecular pro- 
ton transfer in malonaldchydc 47 . Ab initio simulated annealing can be used to 
explore the potential energy landscape and to locate various minima, such as for 
instance done for protonated water clusters 673 . Molecular dynamics simulations of 
the trimcthylaluminum Al(Cfi3)3 have been carried out in order to investigate the 
properties of the gas-phase dimer 29 . The structures and vibrational frequencies 
of tetrathiafulvalcnc in different oxidation states was probed by ab initio molec- 
ular dynamics 324 . Implanted muons in organic molecules (benzene, 3-quinolyl 
nitronyl nitroxide, para-pyridyl nitronyl nitroxidc, phenyl nitronyl nitroxide and 
para-nitrophenyl nitronyl nitroxide) were investigated using approximate ab ini- 
tio path integral simulations that include the strong quantum broadening of the 
muonium 656 > 657 . 

5.9 Chemical Reactions and Transformations 

Early applications of ab initio molecular dynamics were devoted to reactive scat- 
tering in the gas phase such as CH 2 + H 2 -> CH 4 669 or H" + CH 4 -> CH 4 + 
H~ 365 . The "on-the-fly" approach can be compared to classical trajectory cal- 
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culations on very accurate global potential energy surfaces. This was for instance 
done for the well-studied exothermic exchange reaction F + H2 — > HF + H in 
Ref. 565 . Other gas phase reactions studied were Li(2p) + H 2 -» LiH( 1 S) + H( 1 S) 
in Ref. 387 , F + C 2 H 4 -> C 2 H 3 F + H in Ref. 83 , 20 3 -» 30 2 in Ref. 170 , F" + 
CH3CI — » CH 3 F + Cl~ in Ref. 605 , hydroxyl radical with nitrogen dioxide radi- 
cal 165 , formaldehyde radical anion with CH3CI in Ref. 700 , the reduction of OH* 
with 3-hexanone 215 or the hydrolysis (or solvolysis, Sn2 nucleophilic substitution) 
of methyl chloride with water 2 ' 3 . Photoreactions of molecules slowly become ac- 
cessible to ab initio dynamics, such as for instance the cis-trans photoisomerization 
in ethylene 46 , excited-state dynamics in conjugated polymers 71 , bond breaking 
in the Sg ring 562 , transformations of diradicales 195 ' 196 ; or the So — * Si photo 
isomerization of formaldimine 214 . 

In addition to allowing to study complex gas phase chemistry, ab initio molecular 
dynamics opened the way to simulate reactions in solution at finite temperatures. 
This allows liquid state chemistry to take place in the virtual laboratory where 
thermal fluctuations and solvation effects are included. Some applications out of 
this emerging field are the cationic polymerization of 1,2,5-trioxane 146 > 147 ; the 
initial steps of the dissociation of HC1 in water 353 > 354 ; the formation of sulfuric 
acid by letting SO3 react in liquid water or the acid-catalyzed addition of water 
to formaldehyde 422 . 

Proton transfer is a process of broad interest and implications in many fields. 
Intramolecular proton transfer was studied recently in malonaldchydc 695 > 47 ; a Man- 
nich base 182 , and formic acid dimers 427 . Pioneering ab initio molecular dynamics 
simulations of proton and hydroxyl diffusion in liquid water were reported in the 
mid nineties 640 > 641 > 642 . Related to this problem is the auto-dissociation of pure 
water at ambient conditions 628 > 629 . Since recently it became possible to study 
proton motion including nuclear quantum effects 645 i 410 < 412 by using the ab initio 
path integral technique 395,399,644,404 ^ gee gect 4 4 

Ab initio molecular dynamics also allows chemical reactions to take place in solid 
phases, in particular if a constant pressure methodology is used 56 , see Sect. 4.2. 
For instance solid state reactions such as pressure-induced transformations of 
methane 13 and carbon monoxide 54 or the polymerization 57 and amorphization 56 
of acetylene were investigated. 



5.10 Catalysis and Zeolites 

The polymerization of defines is an important class of chemical reactions that is 
operated on the industrial scale. In the light of such applications the detailed un- 
derstanding of these reactions might lead to the design of novel catalysts. Driven 
by such stimulations several catalysts were investigated in detail such as metal 
alkyles 609 , platinum-phospine complexes , or Grubbs' ruthenium -phosphinc 
complexes l , metallocenes 696 . In addition, elementary steps of various chemi- 
cal processes were the focus of ab initio molecular dynamics simulations. Among 
those are chain branching and termination steps in polymerizations 696 , ethylene 
metathesis 1 , "living polymerization" of isoprene with ethyl lithium 522 , Ziegler- 
Natta heterogenous polymerization of ethylene 79 ' 80 ; Reppe carbonylation of Ni 
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CH=CH 2 using C1(C0 ) 2 20 , or Sakakura-Tanaka functionalization 382 . As in the 
real laboratory, side reactions can occur also in the virtual laboratory, such as e.g. 
the /3-hydrogen elimination as an unpredicted reaction path 383 . A digression on 
using finite-temperature ab initio dynamics in homogeneous catalysis research can 
be found in Ref. 697 . 

Zeolites often serve as catalysists as well and are at the same time ideal can- 
didates for finite-temperature ab initio simulations in view of their chemical com- 
plexity A host Of different Studies 559,100,268,614,545,206,560,598,207,315,208,209,546 CQn _ 

tributcd greatly to the understanding of these materials and the processes occurring 
therein such as the initial stages of the methanol to gasoline conversion 599 . Het- 
erogenous catalysts are often poisoned, which was for instance studied in the case 
of hydrogen dissociation on the Pd(100) surface in the presence adsorbed sulfur 
layers 257 . 



5.11 Biophysics and Biochemistry 

Applications of ab initio molecular dynamics to molecules and processes of interest 
in life sciences begin to emerge 18 > 113 . Investigations related to these interests are 
investigations of the crystal structure of a hydrated RNA duplex (sodium guanylyl- 
3'-5'-cytidine nona-hydrate) 311 , structure models for the cytochrom P450 enzyme 
family 547 ' 548 ' 549 , nanotubular polypeptides 112 , a synthetic biomimetic model of 
galactose oxidase 523 , aspects of the process of vision in form of the 11- cis to 
all-trans isomerization in rhodopsin 67,68,474^ m terconversion pathways of the pro- 
tonated /3-ionone Schiff base 615 , or of the binding properties of small molecules 
of physiological relevance such as O2, CO or NO to iron-porphyrines and its com- 
plexes 527,528,529^ 

Proton transport throught water wires is an important biophysical process in 
the chemiosmotic theory for biochemical ATP production. Using the ab initio 
path integral technique 395 , 399, 644,404 the properties of linear water wires with an 
excess proton were studied at room temperature 419 . Amino acids are important 
ingredients as they are the building blocks of polypeptides, which in turn form 
channels and pores for ion exchange. Motivated by their ubiquity, glycine and 
alanine as well as some of their oligopeptides and helical (periodic) polypeptides 
were studied in great detail 323 . 



5.12 Outlook 

Ab initio molecular dynamics is by now not only a standard tool in academic re- 
search but also becomes increasingly attractive to industrial researchers. Analysis 
of data bases, see caption of Fig. 1 for details, uncovers that quite a few companies 
seem to be interested in this methodology. Researchers affiliated to Bayer, Corning, 
DSM, Dupont, Exxon, Ford, Hitachi, Hoechst, Kodak, NEC, Philips, Pirelli, Shell, 
Toyota, Xerox and others cite the Car-Parrinello paper Ref. 108 or use ab initio 
molecular dynamics in their work. This trend will certainly be enhanced by the 
availability of efficient and general ab initio molecular dynamics packages which are 
commercially available. 
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